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Preface 


This document provides a comprehensive description of LSODE, a solver for 
initial value problems in ordinary differential equation systems. It is intended to 
bring together numerous materials documenting various aspects of LSODE, 
including technical reports on the methods used, published papers on LSODE, 
usage documentation contained within the LSODE source, and unpublished notes 
on algorithmic details. 

The three central chapters — on methods, code description, and code usage — are 
largely independent. Thus, for example, we intend that readers who are familiar 
with the solution methods and interested in how they are implemented in LSODE 
can read the Introduction and then chapter 3, Description of Code, without 
reading chapter 2, Description and Implementation of Methods. Similarly, those 
interested solely in how to use the code need read only the Introduction and then 
chapter 4, Description of Code Usage. In this case chapter 5, Example Problem, 
which illustrates code usage by means of a simple, stiff chemical kinetics problem, 
supplements chapter 4 and may be of further assistance. 

Although this document is intended mainly for users of LSODE, it can be used 
as supplementary reading material for graduate and advanced undergraduate 
courses on numerical methods. Engineers and scientists who use numerical 
solution methods for ordinary differential equations may also benefit from this 
document. 
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Chapter 1 

Introduction 


This report describes a FORTRAN subroutine package, LSODE, the Livermore 
Solver for Ordinary Differential Equations, written by Hindmarsh (refs. 1 and 2), 
and the methods included therein for the numerical solution of the initial value 
problem for a system of first-order ordinary differential equations (ODE’s). Such 
a problem can be written as 


i 55 % = 

y(5o) = y 0 = Given * 


(i.i) 


where y, y y, and f are column vectors with N (> 1) components and is the 
independent variable, for example, time or distance. In component form equa- 
tion (1.1) may be written as 


^ = 4*<5> >»(5U) 

y.-(^o) = y>. o = Given 


i = 1 (1.2) 


The initial value problem is to find the solution function y at one or more values 
of ^ in a prescribed integration interval [^o,^ n dL where the initial value of y , y Q 
at ^ is given. The endpoint, may not be known in advance as, for 
example, when asymptotic values of y as £, — » <» are required. 

Initial value, first-order ODE’s arise in many fields, such as chemical kinetics, 
biology, electric network analysis, and control theory. It is assumed that the 
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problem is well posed and possesses a solution that is unique in the interval of 
interest. Solution existence and uniqueness are guaranteed if, in the region of 
interest, f is defined and continuous and for any two vectors y and y* in that 
region there exists a positive constant ££ such that (refs. 3 and 4) 

|f(y.£) - f(y*.4| s se||y - y*|, (i.3) 

which is known as a Lipschitz condition. Here ||*|| denotes a vector norm (e.g., 
ref. 5), and the constant ^ is known as a Lipschitz constant of f with respect to y . 

The right-hand side f of the ODE system must be a function of y and £, only. It 
cannot therefore involve y at previous ^ values, as in delay or retarded ODE’s or 
integrodifferential equations. It cannot also involve random variables, as in 
stochastic differential equations. A second- or higher-order ODE system must be 
reduced to a first-order ODE system. 

The solution methods included in LSODE replace the ODE’s with difference 
equations and then solve them step by step. Starting with the initial conditions at 
So approximations (= Yi.w * = L- JV) to the exact solution y(^„) [=>?,■(£;„), 
i = 1,...,N] of the ODE’s are generated at the discrete mesh points”^ ( n = 1,2,...), 
which are themselves determined by the package. The spacing between any two 
mesh points is called the step size or step length and is denoted by h n , where 

K = ~ Z„- V d-4) 

An important feature of LSODE is its capability of solving “stiff’ ODE problems. 
For reasons discussed by Shampine (ref. 6) stiffness does not have a simple 
definition involving only the mathematical problem, equation (1.1). However, 
Shampine and Gear (ref. 7) discuss some fundamental issues related to stiffness 
and how it arises. An approximate description of a stiff ODE system is that it 
contains both very rapidly and very slowly decaying terms. Also, a characteristic 
of such a system is that the NxN Jacobian matrix J (= 3£/3y ), with element Jy 
defined as 

Jjj = i,j = 1 N, (1.5) 


has eigenvalues {A,,} with real parts that are predominantly negative and also vary 
widely in magnitude. Now the solution varies locally as a linear combination of 
the exponentials {e^ Re ^}, which all decay if all Re(X, ) < 0, where Re(X,) is the 
real part of X,-. Hence for sufficiently large ^ (> l/max|Re(X/)|, where the bars ]*| 
denote absolute value), the terms with the largest Re(X;) will have decayed to 
insignificantly small levels while others are still active, and the problem would be 
classified as stiff. If, on the other hand, the integration interval is limited to 
l/max|Re(X/)|, the problem would not be considered stiff. 
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In this discussion we have assumed that all eigenvalues have negative real 
parts. Some of the Re(^) may be nonnegative, so that some solution components 
are nondecaying. However, the problem is still considered stiff if no eigenvalue 
has a real part that is both positive and large in magnitude and at least one 
eigenvalue has a real part that is both negative and large in magnitude (ref. 7). 
Because the {A,,} are, in general, not constant, the property of stiffness is local in 
that a problem may be stiff in some intervals and not in others. It is also relative in 
the sense that one problem may be more stiff than another. A quantitative 
measure of stiffness is usually given by the stiffness ratio max[-Re(A,/)]/min 
[-Re(^)]. This measure is also local for the reason given previously. Another 
standard measure for stiffness is the quantity max[-Re(A, l )]|2; en <| - ^q|. This 
measure is more relevant than the previous one when |^ en( j - £ol is a better 
indicator of the average “resolution scale” for the problem than l/min[— Re(A,,-)]. 
(In some cases min[-Re(?t f )] = 0.) 

The difficulty with stiff problems is the prohibitive amounts of computer time 
required for their solution by classical ODE solution methods, such as the popular 
explicit Runge-Kutta and Adams methods. The reason is the excessively small 
step sizes that these methods must use to satisfy stability requirements. Because 
of the approximale nature of the solutions generated by numerical integration 
methods, errors are inevitably introduced at every step. For a numerical method 
to be stable, errors introduced at any one step should not grow unbounded as the 
calculation proceeds. To maintain numerical stability, classical ODE solution 
methods must use small step sizes of order l/max[-Re(^)] even after the rapidly 
decaying components have decreased to negligible levels. Examples of the step 
size pattern used by an explicit Runge-Kutta method in solving stiff ODE problems 
arising in combustion chemistry are given in references 8 and 9. Now, the size of 
the integration interval for the evolution of the slowly varying components is of 
order l/min[-Re(A,/)]. Consequently, the number of steps required by classical 
methods to solve the problem is of order max[-Re(^)]/min[-Re(k;)], which is 
very large for stiff ODE’s. 

For stiff problems the LSODE package uses the backward differentiation 
formula (BDF) method (e.g., ref. 10), which is among the most popular currently 
used for such problems (ref. 11). The BDF method possesses the property of stiff 
stability (ref. 10) and therefore does not suffer from the stability step size constraint 
once the rapid components have decayed to negligible levels. Throughout the 
integration the step size is limited only by accuracy requirements imposed on the 
numerical solution. Accuracy of a numerical method refers to the magnitude of 
the error introduced in a single step or, more precisely, the local truncation or 
discretization error. The local truncation error dw at is the difference between 
the computed approximation and the exact solution, with both starting the 
integration at the previous mesh point and using the exact solution y(^ n _i) 
as the initial value. The local truncation error on any step is therefore the error 
incurred on that step under the assumption of no past errors (e.g., ref. 12). 

The accuracy of a numerical method is usually measured by its order. A 
method is said to be of order ^ if the local truncation error varies as /$ +1 . More 
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1. Introduction 

precisely, a numerical method is of order q if there are quantities £ and h 0 (> 0) 
such that (refs. 3 and 13) 

|dj < Ch%+ 1 for all 0 < h n < h Qj (1.6) 

where |dn| is an N-dimensional column vector containing the absolute values of 
the di n 0 = The coefficient vector £ may depend on the function defining 

the ODE and the total integration interval, but it should be independent of the step 
size h n (ref. 13). Accuracy of a numerical method refers to the smallness of the 
error introduced in a single step; stability refers to whether or not this error grows 
in subsequent steps (ref. 7). 

To satisfy accuracy requirements, the BDF method may have to use small step 
sizes of order l/max(Re | Xj\) in regions where the most rapid exponentials are 
active. However, outside these regions, which are usually small relative to the 
total integration interval, larger step sizes may be used. 

The LSODE package also includes the implicit Adams method (e.g., refs. 4 and 
10), which is well suited for nonstiff problems. Both integration methods belong 
to the family of linear multistep methods. As implemented in LSODE these 
methods allow both the step size and the method order to vary (from 1 to 12 for 
the Adams method and from 1 to 5 for the BDF method) throughout the problem. 
The capability of dynamically varying the step size and the method order is very 
important to the efficient use of linear multistep methods (ref. 14). 

The LSODE package consists of 21 subprograms and a BLOCK DATA module. 
The package has been designed to be used as a single unit, and in normal 
circumstances the user needs to communicate with only a single subprogram, also 
called LSODE for convenience. LSODE is based on, and in many ways resembles, 
the package GEAR (ref. 15), which, in turn, is based on the code DIFSUB, written 
by Gear (refs. 10 and 16). All three codes use integration methods that are based 
on a constant step size but are implemented in a manner that allows for the step 
size to be dynamically varied throughout the problem. There are, however, many 
differences between GEAR and LSODE, with the following important 
improvements in LSODE over GEAR; (1) its user interface is much more 
flexible; (2) it is more extensively modularized; and (3) it uses dynamic storage 
allocation, different linear algebra modules, and a wider range of error types (ref. 
17). Most significantly, LSODE has been designed to virtually eliminate the need 
for user adjustments or modifications to the package before it can be used 
effectively. For example, the use of dynamic storage allocation means that the 
required total storage is specified once in the user-supplied subprogram that 
communicates with LSODE; there is no need to adjust any dimension declarations 
in the package. This feature, besides making the code easy to use, minimizes the 
total storage requirements; only the storage required for the user’s problem needs 
to be allocated and not that called for by a code using default values for parameters, 
such as the total number of ODE’s, for example. The many different capabilities 
of the code can be exploited quite simply by setting values for appropriate 
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parameters in the user’s subprogram. Not requiring any adjustments to the code 
eliminates the user’s need to become familiar with the inner workings of the code, 
which can therefore be used as a “black box,” and, more importantly, eliminates 
the possibility of errors being introduced into the modified version. 

The remainder of this report is organized as follows: In chapter 2 we describe 
the numerical integration methods used in LSODE and how they are implemented 
in practice. The material presented in this chapter is based on, and closely 
follows, the developments by Gear (refs. 10 and 18 to 20) and Hindmarsh (refs. 1, 
2, 15, 21, and 22). Chapter 3 describes the features and layout of the LSODE 
package. In chapter 4 we provide a detailed guide to its usage, including possible 
user modifications. The use of the code is illustrated by means of a simple test 
problem in chapter 5. We conclude this report with a brief discussion on code 
availability in chapter 6. 
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Chapter 2 

Description and Implementa- 
tion of Methods 

2.1 Linear Multistep Methods 

The numerical methods included in the packaged code LSODE generate 
approximate solutions X, to the ordinary differential equations (ODE’s) at discrete 
points (n = 1,2,...). Assuming that the approximate solutions Xi -j have been 
computed at the mesh points 2^,_y (j = 1,2,...), these methods advance the solution 
to the current value 2^, of the independent variable by using linear multistep 
formulas of the type 


*1 

Y = Y a.Y n , +h 

— n / j j — n-j n 



j = i 


7=0 


( 2 . 1 ) 


where the current approximate solution vector Xi consists of N components. 





( 2 . 2 ) 


and the superscript T indicates transpose. In equation (2.1), £j_y [= IQLn-j)] is the 
approximation to the exact derivative vector at 2;„_y, y (2^_y) [= f( y (£ n _y))], where 
for notational convenience the 2; argument of f has been droppedfthe coefficients 
{ay} and (Py } and the integers K\ and K 2 are associated with a particular method; 
and h n (= 2; n — 2^,_i) is the step size to be attempted on the current step [2^_i,2^]. 
The method is called linear because the {Xy} and {£} occur linearly. It is called 
multistep because it uses information from several previous mesh points. The 
number max(Jq, K 2 ) gives the number of previous values involved. 

The values /ifi = 1 and = q - 1 produce the popular implicit Adams, or 
Adams-Moulton (AM), method of order q ; 
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<7-1 

y„=y„_ 1 + a„XM«-;- (23) 

;-o 

The method is called implicit because it uses the as yet unknown fj, to compute 
Y n . The method order q means that if equation (2.3) is solved with all past values 
being exact, the resulting Xn will differ from the exact solution y (%*) to the ODE 
system by a local truncation error that is of order 0(// <7+I ) for small values of H- 
max|/i*|. 

The choice K\ = q y K 2 = 0 results in the backward differentiation formula 
(BDF) method of order q: 

Y„=i> ; Y n _,+A n P 0 tr < 2 ' 4 ) 

The term “backward differentiation formula” is used to describe the method 
because equation (2.4), upon division by h n $ o and rearrangement of terms, can be 
regarded as an approximation for y(^«) in terms of Xi> (refs. 15 

and 17). 

The two methods can be written in the general form 

Y„ = + KhU = ln + A„P 0 f(Y n ). (2-5) 

where contains previously computed information and is given by 






;= i 


for the AM method of order q , and 

q 

Yn = X a j^n-j 

j = 1 


(2.6a) 


(2.6b) 


for the BDF method of order q. 

The coefficients {a,} and (P;} are determined such that equations (2.3) and 
(2.4) will be exact if the solution to equation (1.1) is a polynomial of degree q or 
less. Stability characteristics limit q in equation (2.4) to 6 (ref. 10). In LSODE, 
however, BDF’s of order up to only 5 are used because of additional stability 
considerations (refs. 7 and 23). The coefficients (oty) and ( Py 1 for the two 
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methods are given by Gear (ref. 10) for q < 6. In equation (2.5), although the 
subscript n has been attached to the step size h, indicating that h n is the step size to 
be attempted on the current step, the methods used in LSODE are based on a 
constant h. When the step size is changed, the data at the new spacing required to 
continue the integration are obtained by interpolating from the data at the original 
spacing. Solution methods and codes that are based on variable step size have 
also been developed (refs. 17, 23, and 24) but are not considered in the present 
work. 


2.2 Corrector Iteration Methods 

If Po = o, the methods are called explicit because they involve only the known 
values {Y and {£,_;}, and equation (2.1) is easy to solve. If, however, 
po * 0, the methods are called implicit and, in general, solution of equation (2.1) is 
expensive. For both methods, equations (2.3) and (2.4), Po is positive for each q 
and because f is, in general, nonlinear, some type of iterative procedure is needed 
to solve equation (2.5). Nevertheless, implicit methods are preferred because they 
are more stable, and hence can use much larger step sizes, than explicit methods 
and are also more accurate for the same order and step size (refs. 4, 10, and 12). 
Explicit methods are used as predictors, which generate an initial guess for Y„- 
The implicit method corrects the initial guess iteratively and provides a reasonable 
approximation to the solution of equation (2.5). 

The predictor-corrector process for advancing the numerical solution to 
therefore consists of first generating a predicted value, denoted by Xi > anci then 
correcting this initial estimate by iterating equation (2.5) to convergence. That is, 
starting with the initial guess 1^ 0] , approximations Xj, ml (m = 1,2,...,A/) are 
generated (by using one of the techniques discussed below) until the magnitude of 
the difference in two successive approximations approaches zero within a specified 
accuracy. The quantity is the approximation obtained on the mth iteration, 
the integer M is the number of iterations required for convergence, and we accept 
Y^l as an approximation to the exact solution y at and therefore denote it by 
although, in general, it does not satisfy equation (2.5) exactly. 

At each iteration m the quantity /i n Xn"^> which is defined here, is computed 

from x! m] by the relation 


i !," 1 - v „ + M ,^” 1 « 7 ) 

Trwl 

Now, as discussed by Hindmarsh (ref. 21) and shown later in this section, if X« 
converges asm->=», the limit, that is, Hm Y„ , must be a solution of 
equation (2.5) and YM converges to f* [= f(Y„)], the approximation to y (<?„)• 
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Hence h n is the mih estimate for h n in and lim A The predicted 

value of hnf,,, denoted by is also obtained from equation (2.7) (by setting 

m - 0). In practice, we terminate the calculation sequence at a finite number M of 
iterations and accept as an approximation to h ^ the quantity h n Y n m A 
which is obtained from Y^ by using equation (2.7). Note that Y n is only an 
approximation to f* because does not, in general, satisfy equation (2.5) 
exactly (see eqs. (2.5) and (2.7)). Moreover, because Y^ is defined to satisfy 
the solution method, in the sense of equation (2.7), it is not necessarily equal to 
£(y[MJ) Therefore Xn^ and do not necessarily satisfy the ODE, equa- 

tion (1.1). Thus, in practice, to advance the solution, the methods use the {Y y ) (e.g., 
see eqs. (2.8a) and (2.8b)), rather than the {£} as written in equation (2.1). 

After convergence of the estimates xS m] > we could define to be equal to 
f(Y^), so that Xl^ and Y_\!^ satisfy the ODE exactly. However, besides being 
more expensive because it will require one derivative evaluation, performing this 
operation is actually less stable for stiff equations than using equation (2.7) 
(ref. 25). 

The predicted value at X}?\ is generated by a <?th-order explicit formula 
similar to equations (2.3) and (2.4) (refs. 18 and 20): 

xi 01 = Xn-l + K 5>;x n _ ■ (2.8a) 

;=i 

for the AM method of order q and 

xL 0] = i>;x„-, + A n P,*x„_, (2.8b) 

y=i 

for the BDF method of order q. In these two equations X n -j ]S the approximation 
to computed on the step The coefficients [a]} and {pj} are 

selected such that equation (2.8a) or (2.8b) will be exact if the solution to 
equation (1.1) is a polynomial of degree q or less. 

The predictor step for the two methods can be generalized trivially as 

X? = v;. (2.9) 

where \\f* is given by the right-hand sides of equations (2.8a) and (2.8b), 
respectively, for the AM and BDF methods. 
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To correct the initial estimate given by equation (2.9), that is, to solve 
equation (2.5), LSODE includes a variety of iteration techniques — functional, 
Newton-Raphson, and a variant of Jacobi -Newton. 

2.2.1 Functional Iteration 

To derive the functional iteration technique, also called simple iteration 
(refs. 1 1 and 26) and successive substitution (ref. 27), we rewrite equation (2.5) as 
follows: 


X„ = 0(Y„), (2.10) 

where 

*(Y„) = ¥„ + *„M(X n )- (2- 11 ) 

The (m + l)th estimate, (m = 0,1,...,A/-1), is then obtained from 

equation (2.10) by (e.g., ref. 27) 

Y„ [m+1] = ^Y™) = Vj/„ + • (2.12) 

* [W+l] 

Now equation (2.7) gives the following expression for h n Y l n : 


xr 1 = + » 0 Kir n (2.i3) 

Comparing equations (2.12) and (2.13) gives 

K Y n [m+,] = \f(Y' ml ] (2.14) 

for functional iteration. 

We now define the vector function g(y) by 

g(y) = K f(y) + Z , (2.15) 

HO 

which, upon using equation (2.7), gives 


li 



(2.16) 


2. Description and Implementation of Methods 

g(x|'" 1 ) = A„f(YjM) - h n YlT\ 


By using equation (2. 15) we can rewrite the functional iteration equation (2.12) as 
follows: 

Y‘ m+1] = Y l „ m] + P 0 g(x!, ml j. (2.17) 

Finally the combination of equations (2.14) and (2.16) produces the following 
functional iteration procedure for h n Y n : 


• [m+i] 
h n±n 



(2.18) 


Equation (2.17) is simple to use, but it converges only linearly (ref. 27). In 
addition, for successful convergence the step size may be restricted to very small 
values for stiff problems (refs. 4, 10, 12, 26, and 28), as shown here. By using 
equation (2.14) we can rewrite equation (2.16) as 


g(x \ m] ) = ] - /a(xr n ) , (2.i9) 

for m > 1. Hence, equation (2.17) can be rewritten as 


[m+l] 

—n 


In 


[«] 


K Po 


f^xi" 11 ) - ^x^"* -11 ) ■ 


( 2 . 20 ) 


By using the Lipschitz condition, equation (1.3), we get the following relation 
from equation (2.20): 


[m+l] _ Y [m] 

—n —n 



( 2 . 21 ) 


which shows that the iteration converges, that is, the successive differences 

I Y [m+l] 

—n ~ —n 


12 


decrease, only if 


2.2. Corrector Iteration Methods 


KM •= 


( 2 . 22 ) 


Now stiff problems are characterized by, and often referred to as systems with, 
large Lipschitz constants (e.g., refs. 4, 12, and 26), and so equation (2.22) restricts 
the step size to very small values. Indeed, the restriction imposed by this 
inequality on h n is exactly of the same form as that imposed by stability requirements 
on classical methods, such as the explicit Runge-Kutta method (refs. 4 and 26). 
For this reason, when functional iteration is used, the integration method is 
usually said to be explicit even though it is implicit (ref. 17). 

2.2.2 Newton-Raphson Iteration 

Newton-Raphson (NR) iteration, on the other hand, converges quadratically 
and can use much larger step sizes than functional iteration (refs. 27, 29, and 30). 
Rapid improvement in the accuracy of the estimates is especially important 
because the corrector is iterated to convergence. The reason for iterating to 
convergence is to preserve the stability characteristics of the corrector. If the 
correction process is terminated after a fixed number of iterations, the stability 
characteristics of the corrector are lost (refs. 4 and 12), with disastrous consequences 
for stiff problems. 

To derive the NR iteration procedure, we rewrite equation (2.5) as 

E(X„) = Y„ - ln - A„P 0 f(Y„) = 0, (2.23) 

so that solving equation (2.5) is equivalent to finding the zero of R, The quantity 
R0d m] ) is the residual vector on the mth iteration; that is, it is the amount by 
which fails 1° satisfy equation (2.5). To obtain the (m + l)th estimate, we 
expand equation (2.23) in a Taylor series about the mth estimate, neglect the 
second and higher derivatives, and set R(X^ m+l] ) = 0 because we seek a 
that produces this result (e.g., ref. 27). Performing these operations and then 
rearranging terms give the following relation for the NR iteration technique: 



(2,24) 


where the NxN matrix P is given by 

P = 9R/3Y = I -ftJV- (2.25) 
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2. Description and Implementation of Methods 

In equation (2.25), I is the NxN identity matrix and J is the Jacobian matrix, 
equation (1.5). Comparing equations (2.15) and (2.23) shows that 

R(Y) = -p 0 g(Y), (2.26) 

so that equation (2.24) can be rewritten as follows: 

Xr +n = X l n m] + Po P ' 1 i(^ ml ]. (2-27) 

The NR iteration procedure for h n Y„is derived by subtracting equation (2.7) 
from equation (2.13) and then using equation (2.27). The result is 

h n ± l n m+ ' ] = + p_1 g(x« m1 )' (2-28) 


This iteration will converge provided that the predicted value is sufficiently 
accurate (refs. 4 and 12). The prediction method, equation (2.9), provides a 
sufficiently accurate initial estimate that the average number of iterations per step 
is less than 1.5 (ref. 7). In fact, the predictor is generally as accurate as the 
corrector, which is nonetheless needed for numerical stability. However, much 
computational work is required to form the Jacobian matrix and to perform the 
linear algebra necessary to solve equation (2.27). Now, because the Jacobian does 
not appear explicitly in the ODE’s, equation (1.1), or in the solution method, 
equation (2.5), J need not be very accurate. Therefore, for problems in which the 
analytical Jacobian matrix is difficult or impossible to evaluate, a fairly crude 
approximation such as the finite-difference quotient 


h = 


t + t) z f M 


* y j 


ij = 1 


(2.29) 


is adequate. In equation (2.29), A Y } is a suitable increment for thejth component 
of X- 

Inaccuracies in the iteration matrix may affect the rate of convergence of the 
solution but not the solution if it converges (refs. 4 and 21). Hence this matrix 
need only be accurate enough for the iteration to converge. This beneficial fact 
can be used to reduce the computational work associated with linear algebra, as 
described in chapter 3. 
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2.2. Corrector Iteration Methods 

2.2.3 Jacobi-Newton Iteration 

Jacobi-Newton (JN) iteration (ref. 31), also called Jacobi iteration (ref. 32), is 
obtained from Newton-Raphson iteration by neglecting all off-diagonal elements 
of the Jacobian matrix. Hence for JN iteration 




Jo, i*j 
\ tyfiyj, i=j. 


(2.30) 


This technique is as simple to use as functional iteration because it does not 
require any matrix algebra. Also, it converges faster than functional iteration but, 
in general, not as fast as NR iteration. 

A method closely resembling JN iteration is implemented as a separate method 
option in LSODE. It is like JN iteration in that it uses a diagonal approximation D 
to the Jacobian matrix. However, the diagonal elements D „ are, in general, 
different from and are given by the difference quotient 


D u = 


fid + AY) - f.QQ 


A Yj 




(2.31) 


where the increment vector AY = 0. 1 Po g (Yj^)- If J is actually a diagonal matrix, 

Da = Jh + 0(AF, 2 ), but, in general, D /f effectively lumps together the various 
elements {7^} in row i of J. 


2.2.4 Unified Formulation 

The different iteration methods can be generalized by the recursive relations 

Yf +1] = Yj, m] + P 0 p_l g(^ m] j (2.32) 

and 

= h n Y [ n m] + P- I g^ ml j, (2.33) 

where P depends on the iteration method. For functional iteration P = I, and for 
NR and JN iterations P is given by equation (2.25), where J is the appropriate 
Jacobian matrix, equation (1.5), (2.30), or (2.31). 
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2. Description and Implementation of Methods 

The combination of equations (2.32) and (2.33) gives 

Y„ [ ” +l1 - Y‘"l y'” 1 ), (2.34) 

which shows that if converges as m -> », so does . Equation (2.32) 
shows that if Xj/” 1 converges (to X,) as m °°, g(Xii m| ) -» 0, and therefore we 

see from equations (2.15) and (2.16), respectively, (1) that the converged solution 

satisfies equation (2.5) and (2) that Y^ — > f(X) = fn- 

The predictor-corrector methods can be summarized as follows: 

Predictor: 


Y 10] 
— n 



AO] 


h Y [0] 

n n—n 


- 


Po 


(2.35) 


Corrector: 


= hjjYW) - h n Y™ 

Y^rn+l] = Yj-1 + ftp-'gjYW) ► 

k i [ n m+ ' ] = Kilr ] + p _, i(Y' m] ) 


m = 0, 1 


Af - 1. 


(2.36) 


In 


AM] 


Kin = Kil M] - 


(2.37) 
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2.3 Matrix Formulation 


2.3 Matrix Formulation 


The implementation of linear multistep methods is aided by a matrix formulation 
(ref. 21). This formulation, constructed by Gear (ref. 18), is summarized here. 

To solve for X and h n Y„ by using equations (2.35) to (2.37), we need, and 

therefore must have saved, the L = q+ 1 column vectors Xt-l> ^nX n -2’*"» 

and h n Y n _ q for the AM method of order q , or Xt-l> and h n Y n _\ 

for the BDF method of order q . Hence for the AM method of order q we define 
the NxL history matrix w„_j at by 


W„_, = (Y n _ v h n Y n _ v h n Y n _ 2 ,...,h n Y n - q } (2.38a) 


that is, 


v n-\ - 


Y \ y n-\ K^\ y n~\ KY\ t n-2 ■ ■ 

^2,n-l ^n^2,/j-2 • • 


hn^2,n-q 


\?N,n - 1 h n Y N,n-2 ■ • ■ • h n Y N,n-q j 


(2.39) 


The updated matrix 

W„ - (Y n .A„Y„, A„Y n _ 9+1 ) (2.40a) 

is then constructed at each step % n . The predicted matrix wj^ at is given by 


w£» = (Yf 0 U n Y [ n °U n Y„_ I ,....^Y n _ i?+] | (2.41a) 

For the BDF method of order q these matrices take the form 


W„_, = (Y„_ t . h n Y n _ v Y„_ 2 ,...,Y„_J, 

(2.38b) 

w n = (Y n ,A n Y n ,Y„_ 1 ,....Y n _ 9+1 ), 

(2.40b) 
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2. Description and Implementation of Methods 
and 

wL° ] = (2.41b) 

The matrix formulations for wj^ and w„ are derived as follows: Substituting 
the expression for equation (2.8a) or (2.8b), into that for h n equa- 
tion (2.35), and then using equation (2.6a) or (2.6b) give 




7=1 


'i±: 

Po 


h m Y m : + h Y 


n—n-j 


Po 


*n—n-q 


(2.42a) 


for the AM method of order q and 


UT=X 

7=1 


a 7 ~ 


Po 


Pi • 

—n-j + 

Po 


(2.42b) 


for the BDF method of order q. Equations (2.8a) and (2.42a), or (2.8b) and 
(2.42b), that is, the prediction process, can be rewritten as the matrix equation 

w L 01 = w «-i b > ( 2 - 43 > 

where the LxL matrix B depends on the solution method. For the AM method of 
order q y it is given by 


B 


1 

p; 


0 

Pi - P| 
Po 

P2 ~ P2 
Po 


0 0 
1 0 

0 1 


P* . Eld Eld. 0 0 

P<H Po 

P’ 4 0 0 

P<? Po 


0 0^ 
0 0 


0 0 


(2.44a) 


1 


0 1 
0 0 
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and for the BDF method of order q, 


2.3 Matrix Formulation 



1 0 0 ... 0 
0 0 0 . . . o 
0 1 0 . . . 0 

0 0 1 . . . o 


(2.44b) 


a 


<7-1 


a„ 


* 



1 . 

0 0 0 ... 1 
0 0 0 . . . 0 


J 


The corrector equation, equation (2.36), can be expressed in matrix form as 


[w+t] [m] 

+ P 


where the history matrix on the /nth iteration, is given by 

w]" 11 =(X^U n Yt"’U„Y„_ 1 ,...,/ In Y„_ 9+| ) 


(2.45) 


(2.46a) 


for the AM method and by 
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2. Description and Implementation of Methods 


w W = j^Y; m] ,/ In Y' , " 1 > Y n _ 1 ,...,Y n _ <?+1 } (2.46b) 


for the BDF method, le is the L-dimensional vector 


fe = (p o .1.0 0), 


(2.47) 


and P depends on the iteration technique, as described in section 2.2.4. 
The matrix formulation of the methods can be summarized as follows: 

Predictor: 



w «-l B 


(2.48) 


Corrector: 


[ w 1 U f ( V z. v 

gin =*„fX n ~ h n—n 


w [ m+ l] =w [m] + p-1 




m = 0,1,..., A/ — 1. (2.49) 


w = 


(2.50) 


2.4 Nordsieck’s History Matrix 

Instead of saving information in the form w„_ | , equation (2.38a) or (2.38b), 
Gear (ref. 18) suggested making a linear transformation and storing the matrix 
z„_t given by 


«.-I = w n-tQ- ( 251 > 

where the LxL transformation matrix Q is nonsingular. In particular, Q is chosen 
such that the matrix representation suggested by Nordsieck (ref. 33) is obtained: 


z «-l = 


In-lAtn - ,■-£ Y n _„ 



v (<i) 
In-1 > 

) 


(2.52) 
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that is, the NxL matrix z„_i is given by 


2.4 Nordsieck's History Matrix 


f 

Y hn - 1 K\n- 1 

>2,n-l 


y(<7) 

— T M,n-1 

J<?) 


(2.53) 


|^N,n-l K Y N,n - 1 


yfo) 


y 


In equation (2.53), y%_l is the 7 th derivative of the approximating polynomial for 
Y i n _ j. Because scaled derivatives h£xty-\/j\ are used, Q is independent of the 
step size. However, Q depends on the solution method. The N rows of z„_i are 
numbered from 1 to N, so that the /'th row (i = 1,...//) contains the g + 1 scaled 
derivatives of the fth component, Yj in _], of Xn-t- The q + 1 columns are, however, 
numbered from 0 to q , so that the column number corresponds to the order of the 
scaled derivative stored in that column. Thus theyth column (j = 0,1,.. .,< 7 ), which 
we denote by the vector Zn-lti)* contains the vector hfifin-ilj'- The Nordsieck 
matrix formulation of the method is referred to as the “normal form of the 
method” (ref. 10 ). 

Applying the appropriate transformation matrix Q to the predictor equation, 
equation (2.48), gives 

zj ? 1 = wj, 0I Q = w„_jBQ = z„_iQ-'BQ = z„_,A, (2.54) 

where 



is the predicted NxL Nordsieck history matrix at ^ and 
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2. Description and Implementation of Methods 

A = Q-'BQ (2.56) 

The LxL prediction matrix A provides a gth-order approximation to in terms 
of z„_[ and is therefore the lower-triangular Pascal triangle matrix (ref. 10), with 
element Ay given by 




0, i < j 



ij = 0,1 


(2.57) 


where 



is the binomial coefficient, defined as 



i\ 

Mi-W 


(2.58) 


Hence 


r \ o o o o o cr 

11 0 0 

12 1 

13 3 1 0 ..... 

A ss 


I 2 (q-2Xq-l) (g-2Xg-3Xg-4) 

9 2! 3! 

1 a -\ (<? — — 2X<? — 3) 

q 2 ! 3 ! 

, ?(<MX<?-2) 

* 2! 3! 


(9-1) 

2 ! 


(2.59) 


The principal advantage of using the Nordsieck history matrix is that the matrix 
multiplication implied by equation (2.54) can be carried out solely by repeated 
additions, as shown by Gear (ref. 10). Hence computer multiplications are 
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2.4 Nordsieck's History Matrix 

avoided, resulting in considerable savings of computational effort for large 
problems. Also A need not be stored and z n [0] overwrites z n _ h thereby reducing 
memory requirements. 

Because 



(2.60) 


and An = A,o = 1 for all i, the product zA is computed as follows (refs. 10 and 15): 


For k = 0,1, 1, do : 

I For j = q,q- 1, ...,£ + 1, do : 

z i.j + z ij-v i = h...,N. 


(2.61) 


In this equation the subscripts n and n— 1 have been dropped because the z values 
do not indicate any one value of ^ but represent a continuous replacement process. 
At the start of the calculation procedure given by equation (2.61), z = z n _i ; and at 
the end z = z)P\ The arrow denotes the replacement operator, which means 
overwriting the contents of a computer storage location. For example, 


Zi t 3 Zi t 4 + Zi t 3 


means that z f ,4 is added to z/ 3 and the result replaces the contents of the location 
z, 3 . The total number of additions required in equation (2.61) is Nq(q + l)/2. The 
predictor step is a Taylor series expansion about the previous point and is 
independent of both the integration method and the ODE. 

Another important advantage of using Nordsieck’s formulation is that it makes 
changing step size easy. For example, if at ^ the step size is changed from h n to 
rh n > the new history matrix is obtained from 

z n ^z„C, (2.62) 

where the LxL diagonal matrix C is given by 



cO 


(2.63) 




0 


r q 


J 
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2. Description and Implementation of Methods 

The rescaling can be done by multiplications alone, as follows: 


R = 1 

For j = l,...,#, do : 

\R<r-rR 

The corrector equation corresponding to equation (2.49) is given by 
z [m + l] = w ['«+l]Q = wt' n l Q + p-'g(Y[' n l)feQ = z''" 1 +P _l g(Y[ m l)f, 

where z| m] , the Nordsieck history matrix on the mth iteration, is given by 


(2.64) 


(2.65) 


,[m] = 


\\[m] u yW hn. vM !}n_y[m]W 
| L„ - n n±n > 2 \ ~ n ’ ~ n 


V 


V- 


( 2 . 66 ) 


and 


C = kQ ' (2.67) 

is an /.-dimensional vector 

',)■ (168 > 

For the two solution methods used in LSODE the values of f are derived in 
references 21 and 22 and reproduced in tables 2.1 and 2.2. Methods expressed in 
the form of equations (2.54) and (2.65) are better described as multivalue or L- 
value methods than multistep methods (ref. 10) because it is the number L of 
values saved from step to step that is significant and not the number of steps 
involved. 

The two matrix formulations described here are related by the transformation 
equations (2.51), (2.54), and (2.65) and are therefore said to be equivalent 
(ref. 10). The equivalence means that if the step [^n-i4 n ] is taken by the two 
methods with equivalent past values w„_i and z n -l> ^ ' s > related by equa- 
tion (2.51) through Q, then the resulting solutions w„ and z„ will also be related 
by equation (2.51) through Q, apart from roundoff errors (ref. 21). The 
transformation does not affect the stability properties or the accuracy of the 
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2. Description and Implementation of Methods 

TABLE 2.2.— METHOD COEFFICIENTS FOR BACKWARD 
DIFFERENTIATION FORMULA METHOD IN 
NORMAL FORM OF ORDERS 1 TO 6 


<1 

lo 


*2 

‘3 

*4 

*5 

*6 

1 

1 

1 






2 

2 

3 

3 

3 

1 

3 





3 

6 

11 

11 

11 

6 

11 

1 

11 




4 

24 

50 

50 

50 

35 

50 

10 

50 

1 

50 



5 

120 

274 

274 

274 

225 

274 

85 

274 

15 

274 

1 

274 


6 

720 

1704 

1764 

1764 

1624 

1764 

735 

1764 

175 

1764 

21 

1764 

1 

1764 


method, but roundoff properties and computational effort depend on the 
representation used, as discussed by Gear (ref. 10). 

The first two columns of z n and are identical (see eqs. (2.38a), (2.38b), and 
(2.52)), and so { 0 = Po and f \ = 1. For the same reason the corrector iteration 
procedures for Y„ and h n Y n remain unchanged (see eqs. (2.45), (2.47), and 
(2.65)). However, to facilitate estimation of the local truncation error, a different 
iteration procedure than that given by equation (2.65) is used. To derive the new 
formulation, z^ m+1 ^ is written as 


rm+ 1 ] _ [m+1] 

L n ~ L n 


Am) 


+ z 


[«] 


- z 


[m — 1] 


+ ... + z 


[ 1 ] 



+ z 


[ 0 ] 

n 


or 


[m+1] [0] 

z n =Z n 



(2.69) 


Substituting the difference z£ +1] - zj/ 1 obtained from equation (2.65) into equa- 
tion (2.69) produces 


[m+1] [0] 

Z * =Z n 



JO] 




(2.70) 


where e^ m+1 ^ is defined as 
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2.4 Nordsieck’s History Matrix 



It is clear from this equation that 



(2.71) 


(2.72) 


Equation (2.70) can be used to rewrite g(Xjj m] ), equation (2.16), as follows: 



(2.73) 


because = 1. 

Finally, because only the first two columns of z n enter into the solution of equa- 
tion (2.5), the successive corrections can be accumulated and applied to the 
remaining columns of z n after convergence. Clearly, not updating all columns of 
the Nordsieck history matrix after each iteration results in savings of computational 
effort, especially when a high-order method is used and/or the number of ODE’s 
is large. For additional savings of computer time the history matrix is updated 
only if both (1) the iteration converges and (2) the converged solution satisfies 
accuracy requirements. 

The predictor-corrector formulation utilized in LSODE can be summarized as 
follows: 

Predictor: 

,a1 

(2.74) 


Corrector: 


dxlr’l-M nf’i-twiT -si " 1 


[m+l] [m] p . 




v !'* 11 =y !, 01 + foer * 11 




(2.75) 
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2. Description and Implementation of Methods 




Z 


n 



+&,!• 


(2.76) 


2.5 Local Truncation Error Estimate and Control 1 

The local truncation error is defined to be the amount by which the exact 
solution y(£) to the ODE system fails to satisfy the difference equation of the 
numerical method (refs. 4, 10, 12, and 26). That is, for the linear multistep 
methods, equation (2.1), the local truncation error vector^ at is the residual in 
the difference formula when the approximations {Xj} and {£,*} are replaced by the 
exact solution and its derivative. In LSODE, however, the basic multistep 
formula is normalized by dividing it by 


j 

j-0 



Although the corrector convergence test is performed before the local truncation error 
test (which is done only if the iteration converges), we discuss the accuracy test first 
because the convergence test is based on it. 

2 As discussed in chapter 1, another commonly used definition for the local truncation 
error is that it is the error incurred by the numerical method in advancing the approximate 
solution by a single step assuming exact past values and no roundoff errors (refs. 12, 13, 
and 21). That is, d* is the difference between the numerical approximation X* obtained by 
using exact past values (i.e., { y(£„_j)} and ( y (^_y)}) and the exact solution y (^): 

* 

d„=Y*-y(^). (2.77) 

where, ‘for example, 

X* = X a > ) + ^Po ) (2-78) 

7=1 

for the BDF method of order q . For an explicit method the local truncation error given by 
equation (2.77) and that obtained by using the definition given in the text above (i.e., the 
residual of eq. (2.1)) have the same magnitude. However, for an implicit method the two 
quantities are only approximately proportional to one another (ref. 4), although they agree 
asymptotically in the limit of small step size. 
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2.5 Local IVuncation Error Estimate and Control 
for reasons given by Henrici (ref. 29) and Gear (ref. 10); however, see Lambert 
(ref. 4). For example, the BDF method of order q , equation (2.4), can be 
expressed in this form as 


7= 0V 




Po 


X,-; + *„£„• 


(2.79) 


where ccq = ~ L The local truncation error for this method is then given by 



where consists of N components 



(2.80) 


(2.81) 


If we assume that each y, ( i = 1 possesses derivatives of arbitrarily high 
order, each y/(^ n _/) (/ = j = 1, ...,<?) in equation (2.80) can be expanded in a 

Taylor series about Upon collecting terms the resulting expression for d* can 
be stated compactly as 




k = 0 


(2.82) 


where the {Q} are constants (e.g., ref. 10). A method is said to be of order q if 
Q) = C\ = ... = C q - 0, and C q +\ ^ 0. The local truncation error is then given by 


d„ = q+i4r'y (,+1) (6.)+ o (*r 2 ). aw 


where the terms q+i and C q+ iht l y {q+1 \^n) are, respectively, called the error 
constant and the principal local truncation error (ref. 4). In particular, for the BDF 
method of order q in the normalized form given by equation (2.79) (refs. 22 
and 29) 


c<?+1 “ 7+r 


(2.84a) 


For the implicit Adams method of order q in normalized form (ref. 22) 
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2. Description and Implementation of Methods 

C, + H<ofo + U-<o( 4 

where Hq (<?) and ^oO? + 1 ) are, respectively, the zeroth component of the coefficient 
vectors for the AM method in normalized form of orders q and (q + 1). 

The (q + l)th derivative at ^ y is estimated as follows: As discussed 

in section 2.4, at each step the solution method updates the Nordsieck history 
matrix z n : 


Y h Y — Y 
JLn’ n n±n’ 2 !~ n q\ ~ n 


(2.85) 


For either method of order q the last column of z„, £*(<?), contains the vector 
h%Xji q ^qU which is the approximation to y^ q \i sn )/q\. Now the prediction step 
being a Taylor series method of order q does not alter the last column of z n _j, 
namely the vector h%^%l\!q\. Hence the last column of zj[ 01 , contains the 

vector The difference, z^ q) - zj 0] (<?), is given by 


, , [0), . h n v (</) 

l „(4)-Zn (®) = "Tin 

H ' 


K 

q'- 


v ( 9> _ 
— n-l ~ 



( 2 . 86 ) 


by using the mean value theorem for derivatives. However, equation (2.76) gives 
the following expression for ^( 4 ) 

z n W-z l °\q) = <l q e n . (2.87) 


Equating equations (2.86) and (2.87) gives the following approximation for 
h^ +l YJi q+l ^ if higher-order terms are neglected: 


A; +1 ir= ? !f,e„. (2.88) 

Substituting this equation into equation (2.83) and neglecting higher-order terms 
give the following estimate for d*: 

d„ = C ?+ 1 ^e„. (2.89) 

In order to provide for user control of the local truncation error, it is normalized 
by the error weight vector EWT , r , with element EWT / „ defined by 
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2.5 Local Truncation Error Estimate and Control 


EWT- n = RTOL,.(y. n _,|+ ATOL,., (2.90) 

where the user-supplied local relative (RTOL*) and absolute (ATOL;) error toler- 
ances for the ith solution component are discussed in chapter 4. The solution ^ 
is accepted as sufficiently accurate if the following inequality is satisfied: 



^ N 


i 


l 


f 7 


i,n 


EWT, 


< 1. 


i.n J 


(2.91) 


where ||*|| denotes the weighted root-mean-square (rms) norm, which is used for 
reasons discussed by Hindmarsh (ref. 15). Equation (2.91) can be rewritten as 



i N ( 

iy 


e i,n 


N 2 7 


v EWT inj 




by using equation (2.89). If we define the test coefficient x(q,q) as 


*{q,q) 



(2.92) 


(2.93) 


the accuracy test, equation (2.92), becomes 

IM ^ *(<?■<?) 

If we further define the quantity D q by 



the accuracy test reduces to 


D 


? 

< 1. 


(2.94) 


(2.95) 


(2.96) 


The reason for using two variables in the definition for T will become apparent 
when we discuss step size and method order selection in section 2.7. 
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2.6 Corrector Convergence Test and Control 

The test for corrector convergence is independent of both the integration 
method and the iteration technique and is determined by the magnitude of the 
successive differences h n Y [ ™ ] -h n Y [ ™ 11 . To provide for user control of the 
convergence process, the difference h n Y^ 1 -h n is normalized by the 
error weight vector EWT n , equation (2.90). Now, equation (2.33) provides the 

following expression for h n - h n Y^ m ^ : 


h Y W 

ri n± n 


AX n =° n ’ 


(2.97) 


where we have replaced P 1 g(Y£” J ^) by Now, because 


m 

J«+ 1] V sfil 

= 2A 

7=0 

? 

(see eq. (2.71)) and the test on J^l is |e^| < T(^,<?), equation (2.94), the following 
test for convergence 


£ 


m 




i 


i=i 


.[m] 


[ EWT ,'.nJ 


’ T(q,q) 
2(q + 2) 


(2.98) 


is consistent with the local truncation error test. The empirical factor 2 (q + 2) in 
equation (2.98) guarantees that the implicit equation (2.5) is solved to greater 
accuracy than that required of the numerical solution (refs. 22 and 25). 

To increase computational efficiency, especially when the iteration is clearly 
not converging, LSODE uses the following convergence test instead of equa- 
tion (2.98): 


< T(<7 ,<?) 
" 2(9 + 2 ) 


(2.99) 


The quantity E’ m is related to by 


£ , = £ 
m m 


min(l,1.5c^), 


( 2 . 100 ) 
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where 


2.7 Step Size and Method Order Selection and Change 


c' m =max(0-2c m _,,c m ) (2.101) 


and 

1 ( 2 . 102 ) 

is the estimated convergence rate (refs. 22 and 25). Clearly at least two iterations 
are required before c m can be computed. For the first iteration c'm is set equal to 
the last value of c m from the previous step. For the first iteration of the very first 
step and, in the case of NR or JN iteration, after every update of the Jacobian 
matrix, c' m is set equal to 0.7. Equation (2.100) assumes that the iteration 
converges linearly, that is, lim (e m +i/£ m ) = finite constant c, and essentially 

m— 

anticipates the magnitude of E m one iteration in advance (ref. 15). Equation 
(2.101) shows that the convergence rate of the latest iteration is given much more 
weight than that of the previous iteration. The rationale for this decision is 
discussed by Shampine (ref. 25), who examined various practical aspects of 
implementing implicit methods. 


2.7 Step Size and Method Order Selection 
and Change 


Periodically the code attempts to change the step size and/or the method order 
to minimize computational work while maintaining prescribed accuracy. To 
minimize complications associated with method order and step size selection, the 
new order q' is restricted to the values q - 1 >q> and q + 1, where q is the current 
order. For each q f the step size h\q *) that will satisfy exactly the local error bound 
is obtained by assuming that the highest derivative remains constant. The method 
order that produces the largest h f is used on the next step, along with the 
corresponding h\ provided that the h f satisfies certain restrictions described in 
chapter 3. 

For the case q' = q, h\q) is computed by setting Dqih 1 ) (= value of D q for step 
size /z) = 1 (see eq. (2.96)), so that the local accuracy requirement is satisfied 
exactly. Then because d* varies as h # +1 (see eq. (2.83)), we get 


or 




h’(q) 


\ D <U 


<7+1 


(2.103) 
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2. Description and Implementation of Methods 

where r is the ratio of the step size to be attempted on the next step to its current 
value. The subscript “same” indicates that the same order used on the current step 
is to be attempted on the next step. 

For the case q' = q - 1 , d n {q - 1 ) is of order q, where the variable q - 1 indicates 
the method order for which the local truncation error is to be estimated, and 

d n (q-D = C q hly (<, \^ n ), (2.104) 


where C q = | do (q) - do(q - 1)1 for the AM method and \lq for the BDF method 
(refs. 22 and 29). Now, the last column of z„, z n (q), contains the vector h^ q) /ql 
(see eq. (2.85)), and so d n (q - 1) is easily calculated. On using the rms norm, 
equation (2.91), the error test for q' = q - 1 becomes 


i N 


(?) 




cm: 


EWT; 


i,n 


? 

< 1. 


(2.105) 


If we define the test coefficient x(q,q - 1) as MC q q\, equation (2.105) can be 
written as 





x($,9-i) 1) 


? 

< i. 


(2.106) 


where n A (q) is the ith element of z n (q). The first variable in the definition for T 
gives the method order used on the current step. The second variable indicates the 
method order for which the local truncation error is to be estimated. 

The step size h\q- 1 ) to be attempted on the next step, if the order is reduced to 
q — 1, is obtained by using exactly the same procedure that was utilized for the 
case q’ = q , that is, by setting D </ _j(/z r ) = 1. Because &n(q - 1) varies as the 
resulting step size ratio r^own is given by 


h r (q~\) 


'down 


D 


q-\ 


(2.107) 


The subscript “down” indicates that the order is to be reduced by 1 . 
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2,7 Step Size and Method Order Selection and Change 
For the case q' = q + 1 the local truncation error ^{q + 1) is of order q + 2 and is 
given by 


a,(9 + l) = C„ 2 /,f 2 y <, * 2) ft,). (2.108) 

where C q + 2 = | to(q + 2) - (o(q + 1)| for the AM method and 1 l(q + 2) for the BDF 
method (refs. 22 and 29). This case is more difficult than the previous two cases 
because equation (2.108) involves the derivative of order q + 2. The derivative 
y (9+2) (U is estimated as follows. Equation (2.88) shows that the vector Ufa is 
approximately proportional to /t# +I Y/^ +1 Vg!. We difference the quantity 
over the last two steps and use the mean value theorem for derivatives to get 


« 9 Ve„ = H q e n 
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(2.109) 


Hence the error test for q' = q + 1 becomes 


iff 

'c q +2^<?e.n\ 


EWT- 

1 1=1 \ 

, ,n y 


( 2 . 110 ) 


where we have again used the rms norm and is the ith component of Vg n . If 
we define the test coefficient x (q y q + 1) as 1 l(C q + 2 q^q)> the error test, equa- 
tion (2.1 10), can be rewritten as 
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(2.111) 


To solve for h\q + 1), we use the same procedure as for h\q) and h\q - 1). The 
resulting ratio r up is given by 
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1 

q+2 


h^q + 1 ) 

“P h 

n 




( 2 . 112 ) 


The subscript “up” indicates that the order is to be increased by 1. 

After a suitable value for the step size ratio r has been computed, the step size 
h' to be attempted next is calculated: 


h’=rh n . 


(2.113) 


If the step size and/or the method order is changed, the Nordsieck history 
matrix has to be modified. For the case q'=q and A V h n , the q + 1 columns of z n 
are scaled, as described in section 2.4 (see eqs. (2.62) to (2.64)). For the case 
q 1 - q - 1 and h f *h n , the same scaling is performed on the first q columns; the last 
column of the old z n is ignored because it is not needed on subsequent steps. 

If q' = q + 1, z„ must be augmented by a column containing the vector 
(hyt* V(<y + 1)!. The column addition is done in two stages. First, by using 
equation (2.88) we derive the following expression for h% +l Y}p + ty(q + 1)!: 


Y (g+n 

(<7 + 1)! (9 + 1)! q + 1 

and the new column, in{q + 1), is given by 


(2.114) 


^ e 

Z«(9 + 1) = - £z 7- (2-115) 

<7 + 1 

Second, in order to account for any change in the step size, all q + 2 columns of z n 
are rescaled as before. 

Another factor that must be considered if the step size and/or method order is 
changed is that the iteration matrix P, equation (2.25), may be altered even if the 
Jacobian matrix is current. To minimize convergence failures caused by an 
inaccurate P, it must be updated if the coefficient hfo has changed significantly 
since the last evaluation of P. 


2.8 Interpolation at Output Stations 

It frequently happens that the user requires the solution at values ^ Q ut,b 
^out,2>— of the independent variable other than the internally generated {^ n }. It is 
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2.8 Interpolation at Output Stations 
therefore important that in implementing the solution method provision be made 
for the efficient computation of the solution at the required output stations. 
Moreover, the procedure used for these computations should not adversely affect 
the efficiency of the integration beyond the output station. Such a situation arises, 
for example, if the method has to adjust the step size to “hit” the output station 
exactly. Because the Nordsieck history array is used to store past history 
information, the solution can be generated at the output stations quite easily, as 
described next. 

For each ^ out the integration is continued until the first mesh point n for which 
^ and then the solution at i^ u , is obtained by interpolation. Now the 
solution and its scaled derivatives up to order q' n +\ are available at Here q' n +\ 
is the order to be attempted on the next step, that is, [£„-£n+l]- Hence the solution 
at £ out , KUt). is computed by using a (<?;,+ |)th-order Taylor series expansion 
about \ n and is given by 


XM- Y. +(U -%,% + fc™ Y, 

(u -U 


W n+l 


+ . . . + 


(<7 n+l )• 


X? 


n+l) _ <? ^' (^ out ^ n ) 


*=0 


k\ 


Y®. ( 2 . 116 ) 


If we define the quantity r by 


t = 



( 2 . 1 17 ) 


where h' n +\ is the step size to be attempted on the next step, equation ( 2 . 116 ) can 
be rewritten as 



( 2 . 118 ) 


Now 



y(*) 


is the kih column £*(*) of and so equation ( 2 . 118 ) can be expressed compactly 
as 
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0 n+l 

xfcut )=X rt 2«W- 

fc=o 


(2.119) 


Because the solution is accurate to order q’ n +\ at ^ and a (^ r „ + i)th-order Taylor 
series expansion is used to compute Y(^ ut ), the latter is also accurate to order 

q n +\ ■ 

The solution at ^ ut , equation (2.119), can be evaluated by additions and 
multiplications alone by using Homer’s rule (ref. 13): 


>i(Sout) n+l)’ * 

r _ ^out ~ ^>n 
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r*-* r «+i -* 
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(2.120) 


The Taylor series expansion method can be used to compute the solution 
derivative of any order (up to q' n +\) at ^ out . For example, the pth-derivative at 
£out> is given by 


Y^(£ ) = Y f ^ +(£ - £ )y ( ^ +1) + 

— '^out' —n ^Sout n T 


. + 


(^°ut Y ^’" +l ) 

fc.,-.*) 





(* - M 



Y (*) 
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( 2 . 121 ) 


upon using equation (2.117). Substituting = fc!Zfl(£)/(Jr n+ i)* into equa- 
tion (2.121) produces 
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( 2 . 122 ) 
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2.9 Starting Procedure 

2.9 Starting Procedure 

At the outset of the integration, information is available at only the initial point 
£o- Hence multistep methods cannot be used on the first step. The difficulty at 
the initial point is resolved easily by starting the integration with a single-step, 
first-order method. The Nordsieck history matrix zq at £o is constructed from the 
initial conditions yo and the ODE’s as follows: 

z 0 (0)^Y 0 =y o (2.123) 

and 


ZoU^oXo^oHyo^o). (2.124) 


where ho is the step size to be attempted on the first step. 

As the integration proceeds, the numerical solutions generated at the points 
^ 2 *-. provide the necessary values for using multistep methods. Hence, as the 
numerical solution evolves, the method order and step size can be adjusted to 
their optimal values by using the procedures described in section 2.7. 
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Chapter 3 

Description of Code 


3.1 Integration and Corrector Iteration Methods 


The packaged code LSODE has been designed for the numerical solution of a 
system of first-order ordinary differential equations (ODE’s) given the initial 
values. It includes a variable-step, variable-order Adams-Moulton (AM) method 
(suitable for nonstiff problems) of orders 1 to 12 and a variable-step, variable- 
order backward differentiation formula (BDF) method (suitable for stiff problems) 
of orders 1 to 5. However, the code contains an option whereby for either method 
a smaller maximum method order than the default value can be specified. 

Irrespective of the solution method the code starts the integration with a first- 
order method and, as the integration proceeds, automatically adjusts the method 
order (and the step size) for optimal efficiency while satisfying prescribed accuracy 
requirements. Both integration methods are step-by-step methods. That is, 
starting with the known initial condition y(^o) at ^o, where y is the vector of 
dependent variables, % is the independent variable, and is its initial value, the 
methods generate numerical approximations Xn to the exact solution y (^n) at the 
discrete points \ n (n = 1,2,...) until the end of the integration interval is reached. 
At each step &-1&] both methods employ a predictor-corrector scheme, wherein 
an initial guess for the solution is first obtained and then the guess is improved 
upon by iteration. That is, starting with an initial guess, denoted by Xn . 
successively improved estimates XT ' (m = are generated until the iteration 

converges, that is, further iteration produces little or no change in the solution. 
Here is the approximation computed on the mth iteration, and M is the 
number of iterations required for convergence. 

A standard explicit predictor formula— a Taylor series expansion method devised 
by Nordsieck (ref. 33)— is used to generate the initial estimate for the solution. A 
range of iteration techniques for correcting this estimate is included in LSODE. 
Both the basic integration method and the corrector iteration procedure are identified 
by means of the method flag MF. By definition, MF has the two decimal digits 


METH and MITER, and 

PRfeCStMNS PAGE 


BLANK not filmed 

PAGE jKL- INTENTIONALLY BLANK 
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3. Description of Code 


TABLE 3.1. —SUMMARY OF INTEGRATION METHODS INCLUDED IN LSODE 
AND CORRESPONDING VALUES OF METH, 

THE FIRST DECIMAL DIGIT OF MF 


METH 

Integration method 

1 

2 

Variable-step, variable-order, implicit Adams method of orders I to 12 
Variable-step, variable-order, implicit backward differentiation formula 
method of orders 1 to 5 


TABLE 3.2.— CORRECTOR ITERATION TECHNIQUES AVAILABLE IN LSODE 
AND CORRESPONDING VALUES OF MITER, 

THE SECOND DECIMAL DIGIT OF MF 


MITER 

Corrector iteration technique 

0 

Functional iteration 

I 

Modified Newton iteration with user-supplied analytical Jacobian 

2 

Modified Newton iteration with internally generated numerical Jacobian 

3 

Modified Jacobi-Newton iteration with internally generated numerical 
Jacobian 3 

b 4 

Modified Newton iteration with user-supplied banded Jacobian 

b 5 

Modified Newton iteration with internally generated banded Jacobian 


a Modified Jacobi-Newton iteration with user-supplied analytical Jacobian can be 
performed by specifying MITER = 4 and ML = MU = 0 b (i.e., a banded Jacobian 
with bandwidth of 1). 

^The user must specify the lower (ML) and upper (MU) half-bandwidths of the 
Jacobian matrix. 


MF = 1 0 x METH + MITER, (3.1) 

where the integers METH and MITER indicate, respectively, the integration 
method and the corrector iteration technique to be used on the problem. Table 3. 1 
summarizes the integration methods included in LSODE and the appropriate 
values for METH. The legal values for MITER and their meanings are given in 
table 3.2. The iteration procedures corresponding to MITER = 1 to 5 are 
described as modified Newton iteration techniques because the Jacobian matrix is 
not updated at every iteration. 


3.2 Code Structure 


The double-precision version of the LSODE package consists of the main core 
integration routine, LSODE, the 20 subprograms CFODE, DAXPY, DDOT, 
DGBFA, DGBSL, DGEFA, DGESL, DSCAL, D1MACH, EWSET, IDAMAX, 
INTDY, PREPJ, SOLSY, SRCOM, STODE, VNORM, XERRWV, XSETF, and 
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3.2 Code Structure 


XSETUN, and a BLOCK DATA module for loading some variables. The single- 
precision version contains the main routine, LSODE, and the 20 subprograms 
CFODE, EWSET, INTDY, ISAMAX, PREPJ, R1MACH, SAXPY, SDOT, 
SGBFA, SGBSL, SGEFA, SGESL, SOLSY, SRCOM, SSCAL, STODE, 
VNORM, XERRWV, XSETF, and XSETUN. The subprograms DDOT, 
D1MACH, ID AM AX, ISAMAX, R1MACH, SDOT, and VNORM are function 
routines — all the others are subroutines. The subroutine XERRWV is machine 
dependent. In addition to these routines the following intrinsic and external 
routines are used: DABS, DFLOAT, DMAX1, DMIN1, DSIGN, and DSQRT 
by the double-precision version; ABS, AMAX1, AMIN1, FLOAT, SIGN, and 
SQRT by the single-precision version; and MAXO, MINO, MOD, and WRITE 
by both versions. 

Table 3.3 lists the subprograms in the order that they appear in the code and 
briefly describes each subprogram. Among these, the routines DAXPY, DDOT, 
DGBFA, DGBSL, DGEFA, DGESL, DSCAL, IDAMAX, ISAMAX, SAXPY, 
SDOT, SGBFA, SGBSL, SGEFA, SGESL, and SSCAL were taken from the 
LINPACK collection (ref. 34). The subroutines XERRWV, XSETF, and 
XSETUN, as used in LSODE, constitute a simplified version of the SLATEC 
error-handling package (ref. 35). 

The structure of the LSODE package is illustrated in figure 3.1, wherein a line 
connecting two routines indicates that the lower routine is called by the upper one. 
For subprograms that have different names in the different versions of the code, 
both names are given, with the double-precision version name listed first. Also, 
the names in brackets are dummy procedure names, which are used internally and 
passed in call sequences. The routine Fisa user-supplied subroutine that computes 
the derivatives dy-Jd % (i = 1,...,jV), where y/ is the ith component of y and N is the 
number of ODE’s. Finally, the user-supplied subroutine JAC computes the 
analytical Jacobian matrix J (= di/d y), where f = dy/dt , . 

The code has been arranged as much as possible in a “modular” fashion, with 
different subprograms performing different tasks. Hence the number of 
subprograms is fairly large. However, this feature aids in both understanding and, 
if necessary, modifying the code. To enhance the user’s understanding of the 
code, it contains many comment statements, which are grouped together in blocks 
and describe both the task to be performed next and the procedure to be used. In 
addition, each subprogram includes detailed explanatory notes, which describe 
the function of the subprogram, the means of communication (i.e., call sequence 
and/or common blocks), and the input and output variables. 

Each subprogram contains data type declarations for all variables in the routine. 
Such declarations are useful for debugging and provide a list of all variables that 
occur in a routine. This list is useful in overlay situations. For each data type the 
variables are usually listed in the following order: variables that are passed in the 
call sequence, variables appearing in common blocks, and local variables, in 
either alphabetical order or the order in which they appear in the call sequence and 
the common blocks. 
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TABLE 3.3.— DESCRIPTION OF SUBPROGRAMS USED IN LSODE 


Subprogram 

Description 

Double- 

precision 

version 

Single- 

precision 

version 


LSODE 

LSODE 

Main core integration routine. Checks legality of input, 
sets work array pointers, initializes work arrays, com- 
putes initial integration step size, manages solutions 
of ODE’s, and returns to calling routine with solution 
and errors. 

INTDY 

INTDY 

Computes interpolated values of the specified derivative 
of the dependent variables. 

STODE 

STODE 

Advances the solution of the ODE’s by one integration 
step. Also, computes step size and method order to be 
attempted on the next step. 

CFODE 

CFODE 

Sets method coefficients for the solution and test con- 
stants for local error test and step size and method order 
selection. 

PREPJ 

PREPJ 

Computes the iteration matrix and either manages the 
subprogram call for its LU-decomposition or computes 
its inverse. 

SOLSY 

SOLSY 

Manages solution of linear system arising from chord 
iteration. 

EWSET 

EWSET 

Sets the error weight vector. 

VNORM 

VNORM 

Computes weighted root-mean-square norm of a vector. 

SRCOM 

SRCOM 

Saves and restores contents of common blocks LSOOOI 
and EHOOOl. 

D1MACH 

RIMACH 

Computes unit roundoff of the computer. 

XERRWV 

XERRWV 

Handles error messages. 

XSETF 

XSETF 

Resets print control flag. 

XSETUN 

XSETUN 

Resets logical unit number for error messages. 

DGEFA 

SGEFA 

Performs LU-decomposition of a full matrix by Gaussian 
elimination. 

DGESL 

SGESL 

Solves a linear system of equations using a previously 
LU-decomposed full matrix. 

DGBFA 

SGBFA 

Performs LU-decomposition of a banded matrix by 
Gaussian elimination. 

DGBSL 

SGBSL 

Solves a linear system of equations using a previously 
LU-decomposed banded matrix. 

DAXPY 

SAXPY 

i 

Forms the sum of one vector and another times a 
constant. 

DSCAL 

SSCAL 

Scales a vector by a constant. 

DDOT 

SDOT 

Computes dot product of two vectors. 

IDAMAX 

ISAM AX 

Identifies vector component of maximum absolute value. 
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3. Description of Code 

3.3 Internal Communication 


Communication between different subprograms is accomplished by means of 
both call sequences and the two common blocks EH0001 and LS0001. The 
reason for using common blocks is to avoid lengthy call sequences, which can 
significantly deteriorate the efficiency of the program. However, common blocks 
are not used for variables whose dimensions are not known at compilation time. 
Instead, to both eliminate user adjustments to the code and minimize total storage 
requirements, dynamic dimensioning is used for such variables. 

The common blocks, if any, used by each subprogram are given in tables 3.4 
and 3.5 for the double- and single-precision versions, respectively. These tables 
also list all routines called and referenced (e.g., an external function) by each 
subprogram. Also, to facilitate use of LSODE in overlay situations, all routines 
that call and reference each subprogram are listed. Finally, for each subprogram 
the two tables give dummy procedure names (which are passed in call sequences 
and therefore have to be declared external in each calling and called subprogram) 
in brackets. 

The variables included in the two common blocks and their dimensions, if 
different from unity, are listed in table 3.6. The common blocks contain variables 
that are (1) local to any routine but whose values must be preserved between calls 
to that routine and (2) communicated between routines. The structure of the block 
LS0001 is as follows: All real variables are listed first, then all integer variables. 
Within each group the variables are arranged in the following order: (1) those 
local to subroutine LSODE, (2) those local to subroutine STODE, and (3) those 
used for communication between routines. It must be pointed out that not all 
variables listed for a given common block are needed by each routine that uses it. 
For this reason some subprograms may use dummy names, which are not listed in 
table 3.6. 

To further assist in user understanding and modification of the code, we have 
included in table 3.6 the names of all subprograms that use each common block. 
For the same reason we provide in tables 3.7 and 3.8 complete descriptions of the 
variables in EH0001 and LS0001, respectively. Also given for each variable are 
the default or current value, if any, and the subprogram (or subprograms) where it 
is set or computed. The length LENWM of the array WM in table 3.8 depends on 
the iteration technique and is given in table 3.9 for each legal value of MITER. 

3.4 Special Features 

The remainder of this chapter deals with the special features of the code and its 
built-in options. We also describe the procedure used to advance the solution by 
one step, the corrective actions taken in case of any difficulty, and step size and 
method order selection. In addition, we provide detailed flowcharts to explain the 
computational procedures. We conclude this chapter with a brief discussion of the 
error messages included in the code. 
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TABLE 3.4.— ROUTINES WITH COMMON BLOCKS, SUBPROGRAMS, AND 
CALLING SUBPROGRAMS IN DOUBLE-PRECISION 
VERSION OF LSODE 


Subprogram 
[Dummy 
procedure name] 

Common blocks 
used 

Subprograms 
called and 
referenced 

Calling 

subprograms 

LSODE 

LS0001 

D1MACH EWSET 




F INTDY JAC 




PREPJ SOLSY 




STODE VNORM 




XERRWV 


CFODE 



STODE 

DAXPY 



DGBFA DGBSL 




DGEFA DGESL 

DDOT 



DGBSL DGESL 

DGBFA 


DAXPY DSCAL 

PREPJ 



IDAMAX 


DGBSL 


DAXPY DDOT 

SOLSY 

DGEFA 


DAXPY DSCAL 

PREPJ 



IDAMAX 


DGESL 


DAXPY DDOT 

SOLSY 

DSCAL 



DGBFA DGEFA 

D1MACH 



LSODE 

EWSET 



LSODE 

IDAMAX 



DGBFA DGEFA 

INTDY 

LS0001 

XERRWV 

LSODE 

PREPJ 

LS0001 

DGBFA DGEFA 

STODE 

[PJAC] 


F JAC VNORM 


SOLSY 

LS0001 

DGBSL DGESL 

STODE 

[SLVS] 




SRCOM 

EH0001 LS0001 



STODE 

LS0001 

CFODE F JAC 

LSODE 



PREPJ SOLSY 




VNORM 


VNORM 



LSODE PREPJ 




STODE 

XERRWV 

EH0001 


LSODE INTDY 

XSETF 

EHOOOI 



XSETUN 

EH0001 



BLOCK DATA 

EHOOOI LS0001 
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TABLE 3.5.— ROUTINES WITH COMMON BLOCKS, SUBPROGRAMS, AND 
CALLING SUBPROGRAMS IN SINGLE-PRECISION 
VERSION OF LSODE 


Subprogram 
[Dummy 
procedure name] 

Common blocks 
used 

Subprograms 
called and referenced 

Calling 

subprograms 

LSODE 

LS0001 

EWSET F INTDY 




JAC PREPJ 




RIMACH SOLSY 




STODE VNORM 




XERRWV 


CFODE 



STODE 

EWSET 



LSODE 

INTDY 

LS0001 

XERRWV 

LSODE 

ISAM AX 



SGBFA SGEFA 

PREPJ 

LS0001 

F JAC SGBFA 

STODE 

[PJAC] 


SGEFA VNORM 


RIMACH 



LSODE 

SAXPY 



SGBFA SGBSL 




SGEFA SGESL 

SDOT 



SGBSL SGESL 

SGBFA 


ISAMAX SAXPY 

PREPJ 



SSCAL 


SGBSL 


SAXPY SDOT 

SOLSY 

SGEFA 


ISAMAX SAXPY 

PREPJ 



SSCAL 


SGESL 


SAXPY SDOT 

SOLSY 

SOLSY 

LS0001 

SGBSL SGESL 

STODE 

[SLVS] 




SRCOM 

EH0001 LS0001 



SSCAL 



SGBFA SGEFA 

STODE 

LS0001 

CFODE F JAC 

LSODE 



PREPJ SOLSY 




VNORM 


VNORM 



LSODE PREPJ 




STODE 

XERRWV 

EH0001 


LSODE INTDY 

XSETF 

EH0001 



XSETUN 

EH0001 




48 










TABLE 3.6.— COMMON BLOCKS WITH VARIABLES AND 
SUBPROGRAMS WHERE USED 


Common 

block 

Variables (dimension) 

Subprograms where 
used 

EH0001 

MESFLG LUNIT 

SRCOM XERRWV 
XSETF XSETUN 
BLOCK DATA a 

LS0001 

CONIT CRATE EL(13) 

LSODE INTDY 


ELCCK13, 12) HOLD RMAX 

PREPJ SOLSY 


TESCO(3, 12) CCMAX ELO 

SRCQM STODE 


H HMIN HMXI HU RC TN 
UROUND ILLIN INIT LYH 
LEWT LACOR LSAVF LWM 
LIWM MXSTEP MXHNIL 
NHNIL NTREP NSLAST 
NYH IALTH IPUP LMAX 
MEO NQNYH NSLP ICF 
IERPJ IERSL JCUR JSTART 
KFLAG L METH MITER 
MAXORD MAXCOR MSBP 
MXNCF N NQ NST NFE 
NJE NQU 

BLOCK DATA* 


“Double-precision version only. 


TABLE 3.7 — DESCRIPTION OF VARIABLES IN COMMON BLOCK EH0001, 
THEIR CURRENT VALUES, AND SUBPROGRAMS WHERE THEY ARE SET 


Variable 

Description 

Current 

value 

Subprogram where 
variable is set 

MESFLG 

Integer flag, which controls 
printing of error messages from 
code and has following values 
and meanings: 

0 No error message is printed. 

1 All error messages are printed. 

1 

BLOCK DATA in 
double -precision version 
and XERRWV in single- 
precision version 

LUNIT 

Logical unit number for messages 
from code 

6 

BLOCK DATA in 
double-precision version 
and XERRWV in single- 
precision version 
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TABLE 3.8 — DESCRIPTION OF VARIABLES IN COMMON BLOCK LS0001, THEIR 
CURRENT VALUES, IF ANY, AND SUBPROGRAMS WHERE 
THEY ARE SET OR COMPUTED* 


Variable 

Description 

Current value. 

Subprograms where 



if any 

variable is set or 




computed 

CONIT 

Empirical factor, 0.5/(NQ + 1 ) 

— — 

STODE 


used in convergence test (see 
eq. (2.99)) 



CRATE 

Estimated convergence rate of 

— 

STODE 


iteration 



EL 

Method coefficients in normal 


STODE 


form {^} (see eq. (2.68)), for 
current method order 



ELCO 

Method coefficients in normal 

— 

CFODE 


form for current method of 
orders 1 to MAXORD 



HOLD 

Step size used on last success- 

— 

STODE 


ful step or attempted on last 
unsuccessful step 

Normally 10; 10 4 for very 


RMAX 

Maximum factor by which step 

STODE 


size will be increased when 

first step size increase 



step size change is next 

for problem if no dif- 



considered 

ficulty encountered; 2 
after a failed converg- 
ence or local error test 


TESCO 

Test coefficients for current 


CFODE 


method of orders 1 to 
MAXORD; used for testing 
convergence and local 
accuracy and selecting new 
step size and method order 



CCMAX 

Maximum relative change 

0.3 

LSODE 


allowed in HxELO before 
Jacobian matrix is updated 



ELO 

Method coefficient Gg (see 


STODE 


eq. (2.68)) for current method 
and current order 



H 

Step size either being used on 


LSODE 


this step or to be attempted 
on next step 

0.0 

STODE 

HMIN b 

Minimum absolute value of step 



size to be used on any step 

0.0 


HMXI b 

Inverse of maximum absolute 

LSODE 


value of step size to be used 
on any step 


STODE 

HU 

Step size used on last success- 



ful step 




"Note that some variables appear in the table before they are defined. 

b Default value for this variable can be changed by the user, as described in table 4.6. 
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Variable 


Description 


TABLE 3.8. — Continued. 


Subprograms where 
variable is set or 
computed 


Current value, 
if any 


UROUND 


Relative change in HxELO 
since last update of Jacobian 
matrix 

Value of independent variable 
to which integrator either has 
successfully advanced solution 
or will do so after next step 
Unit roundoff of computer 


Number of consecutive times 
LSODE has been called with 
illegal input for current 
problem 


Integer flag (= 0 or 1) that 
denotes if initialization of 
LSODE has been performed 
(INIT m 1) or not (INTT * 0) 
Base address for Nordsicck 
history array YH of length 
NYHx(MAXORD + I) 

Base address for error weight 
vector EWT of length N 
Base address for array ACOR 
(of length N) containing local 
errors on last successful step 
Base address for an array 
SAVF (of length N), used for 
temporary storage 
Base address for array WM (of 
length LENWM 0 ), required 
for linear algebra associated 
with Jacobian and iteration 


LWM + LENWM 0 


LEWT + 2N 


LEWT + N 


LYH + 

NYHx(MAXORD + 1) 



D1MACH in 
double-precision 
version and 
R1MACH in 
single-precision 
version 
Initialized in 
BLOCK DATA 
(double -precision 
version) and 
LSODE (single- 
precision version). 
Updated in LSODE 
in both versions. 

LSODE 


c The length LENWM of the array WM depends on the iteration technique and is given in 
table 3.9. 
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TABLE 3.8. — Continued. 


Variable 


Description 


Current value, 
if any 


Subprograms where 
variable is set or 
computed 


LIWM 

MXSTEph 

MXHNIL b 


NHNIL 


NTREP 


NSLAST 


NYH 


IALTH 


Base address for integer work 
array IWM 

Maximum number of steps 
allowed on any one call to 
LSODE 

Maximum number of times that 
warning message that step 
size is so small that TNV 
H = TN for next step is 
printed 

Number of times th§t this dif- 
ficulty with small step size 
has been encountered so far 
for problem . 

Number of consecutive times an 
initialization or "first" call 
(see table 4.3) has been made 
to LSODE with same initial 
and final values for integra- 
tion interval 


Number of steps used for 
problem prior to current call 
to LSODE; used to check that 
the limit of MXSTEP steps is 
not exceeded 

Maximum number of ODE's to 
be solved for current problem 
(This number is equal to the 
number of ODE’s specified on 
first call to LSODE.) 

Integer counter, related to step 
size and method order 
changes, with following 
values and meanings: 

0 Select optimal step size and 
method order. 

1 If NQU < MAXORD, save 
vector e (see eqs, (2.76) and 
(2.11 1)) so that an order 
increase can be considered 
on the next step. 

>1 Neither of these two oper- 
ations is to be performed. 


1 

500 

10 


LSODE 

LSODE 

LSODE 

LSODE 


Initialized in 
BLOCK DATA 
(double-precision 
version) and 
LSODE (single- 
precision 
version). Updated 
in LSODE in 
both versions. 

LSODE 


LSODE 


STODE 


b Default value for this variable can be changed by the user, as described in table 4.6. 
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TABLE 3.8. — Continued. 


Variable 


Description 


Current value, 
if any 


Subprograms where 
variable is set or 
computed 


IPUP 


LMAX 

MEO 

NQNYH 

NSLP 

ICF 


IERPJ 


IERSL 


Integer flag, related to Jacobian 
matrix update, with following 
values and meanings: 

0 Jacobian matrix is either 
not needed or does not have 
to be updated. 

>0 Jacobian matrix must be 
updated before corrector 
iteration. 

Maximum number of columns 
of Nordsieck history array 
Integration method specified on 
previous call to LSODE 
Number of dements of 
Nordsieck history array that 
are changed by predictor 
Step number when Jacobian 
matrix was last updated 
An integer flag, related to iter- 
ation convergence, with fol- 
lowing values and meanings: 

0 Solution converged. 

1 Convergence test failed and 
Jacobian matrix is not 
current. 

2 Convergence test failed and 
Jacobian matrix is either 
current or not needed. 

Integer flag, related to singular- 
ity of iteration matrix, with 
following values and 
meanings: 

0 Iteration matrix was suc- 
cessfully LU -decomposed 
(MITER = 1,2, 4, or 5) or 
inverted (MITER = 3) (see 
table 3.2) 

1 Iteration matrix was found 
to be singular. 

Integer flag, related to singular- 
ity of interadon matrix modi- 
fied to account for new 
(HxELO) for MITER = 3 (see 
table 3.2). IERSL has fol- 
lowing values and meanings: 

0 Modified iteration matrix 
was successfully inverted 
and corrections computed. 

1 New matrix was found to be 
singular. 


MAXORD + 1 


NQxNYH 


STODE 


STODE 

STODE 

STODE 

STODE 

STODE 


PREPJ 


SOLSY 
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TABLE 3.8. — Continued. 


Variable 

Description 

Current value. 

Subprograms where 



if any 

variable is set or 




computed 

JCUR 

Integer flag, related to state of 


PREPJ 


Jacobian matrix, with fol- 
lowing values and meanings; 
0 Jacobian matrix is not cur- 


STODE 


rent and may need to be 
updated lata*. 




1 Matrix is current 



JSTART 

Integer flag, used to comm uni- 


LSODE 


cate state of calculation to 
STODE, with following 
values and meanings: 

0 This is the first step for the 
problem. 

! Continue normal calculation 


STODE 


of problem. (This is the value 
returned by STODE to 
facilitate continuation.) 




-1 Take the next step with new 




values for H, MAXORD, N, 
METH (see table 3.1), 
MITER (see table 3.2), 
and/or matrix parameters. 



KFLAG 

A completion code from 
STODE with following values 
and meanings: 

0 Step was successful. 

-1 Requested local accuracy in 


STODE 


solution could not be 
achieved. 




-2 Repeated convergence test 




failures occurred. 



L 

I Number of columns of 

NQ+ 1 

STODE 


Nordsieck array 



METH 

Integration method to be used 


LSODE 


on next step 



MITER 

Iteration technique to be used 


LSODE 

MAXORD 6 

on next step 

Maximum method order to be 

12 for Adams-Moulton 

LSODE 


used for problem 

method and 5 for 
backward different- 
iation formula method 


MAXCOR 

Maximum number of corrector 

3 

LSODE 


iterations to be attempted on 
any one step 



MSBP 

Maximum number of steps for 

20 

LSODE 


which same Jacobian matrix is 
used 




kDefauit value for this variable can be changed by the user, as described in table 4.6. 
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TABLE 3.8. — Concluded. 


Variable 

Description 

Current value, 
if any 

Subprograms where 
variable is set or 
computed 

MXNCF 

Maximum number of corrector 

10 

LSODE 


convergence failures allowed 




on any one step 



N 

Number of ODE’s to be solved 

— 

— 


on next step 



NQ 

Method order either being tried 

— 

STODE 


on this step or to be attempted 




on next step 



NST 

Total number of integration 

— 

LSODE 


steps used so far for problem 


STODE 

NFE 

Total number of derivative 


LSODE 


evaluations required so far for 


STODE 


problem 



NJE 

Total number of Jacobian 

— 

LSODE 


matrix evaluations (and 


PREPJ 


iteration matrix LU- 




decompositions or inversions) 




required so far for problem 



NQU 

Method order used on last suc- 


STODE 


cessful step. 




TABLE 3.9.— LENGTH LENWM 
OF ARRAY WM IN TABLE 3.8 
FOR ITERATION TECHNIQUES 
INCLUDED IN CODE 


MITLR a 

LENWM 1 ’ 

0 

0 

1,2 

N 2 + 1 

3 

N + 2 

4,5 

(2ML + MU + l)N + 2 


a See table 3.2 for description of 
MITER. 

^ is the number of ODE’s and ML 
and MU are defined in table 3.2. 
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3. Description of Code 

The main routine, LSODE, controls the integration and serves as an interface 
between the calling subprogram and the rest of the package. A flowchart of this 
subroutine is given in figure 3.2. In this figure ITASK and ISTATE are user- 
specified integers that specify, respectively, the task to be performed and the state 
of the calculation, that is, if the call to LSODE is the first one for the problem or a 
continuation; if the latter, ISTATE further indicates if the continuation is a normal 
one or if the user has changed one or more parameters since the last call to 
LSODE (see chapter 4 for details). On return from LSODE the value of ISTATE 
indicates if the integration was performed successfully, and if not, the reason for 
failure. The integer JSTART is an internally defined variable used for 
communicating the state of the calculation with the routine STODE. The variables 
T (= £), H, and X are, respectively, the independent variable, the step size to be 
attempted on the next step, and the numerical solution vector. TOUT is the 
value at which the solution is next required. Finally, TCRIT is the £ value that 
the integrator must not overshoot. This option is useful if a singularity exists at or 
beyond TCRIT and is discussed further in chapter 4. 

The subroutine STODE advances the numerical solution to the ODE’s by a 
single integration step It also computes the method order and step size 

to be attempted on the next step. The efficiency of the integration procedure is 
increased by saving the solution history, which is required by the multistep 
methods used in the code, in the form suggested by Nordsieck (ref. 33). The 
Nx(q + 1) Nordsieck history matrix z„_i at \ n -\ contains the numerical solution 
Xi-i and the q scaled derivatives h } n X$-\lj\ (j = l>...,tf), where h n (= ^ n -^n-\) and 
q are, respectively, the current step size and method order and X® = d^Y/d^. 

The flowchart of STODE is presented in Figure 3.3. In this Figure NCF is the 
number of corrector convergence failures on the current step, KFLAG is an 
internally deFined integer used for communication with LSODE, NQ (= q) is the 
method order to be attempted on the current step, and the integer counter IALTH 
indicates how many more steps are to be taken with the current step size and 
method order. The (NQ + l)-dimensional vector { contains the method coefficients 
and depends on both the integration method and the method order; ?o is the zeroth 
component of (! (see eq. (2.68)). The matrix zj°i is the predicted Nordsieck 
history matrix at and the NxN iteration matrix P is given by equation (2.25). 
The variable R is the ratio of the step size to be attempted next to its current 
value, RMAX is the maximum R allowed when a step size change is next 
considered, and HMIN and HMAX are user-supplied minimum and maximum 
absolute values for the step size to be tried on any step. The ratios RHDN, 
RHSM, and RHUP are factors by which the step size can be increased if the new 
method order is NQ - 1, NQ (the current value), and NQ + 1, respectively. 
Finally, NQMAX is the maximum method order that may be attempted on any 
step, and the vector & n (=h n Y n - h is proportional to the local truncation 
error vector at E , n (see eqs. (2.87) an S (2.89)). 
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3.4 Special Features 

3.4.1 Initial Step Size Calculation 

An important feature of LSODE is that it will compute the step size ho to be 
attempted on the first step if the user does not provide a value for it. The 
calculation procedure attempts to produce an ho such that the numerical solution 
Xi generated at the first internal mesh point £,\ will satisfy the local error test. 
Now with either solution technique the code starts the integration with a first- 
order method. Hence the asymptotic local truncation error d K \ in the ith solution 
component at will be equal to (l/2)h\yj(^\) for both the AM and BDF methods 
of order 1 . Here h\ is the step size successfully used on the first step, and is 
the second derivative of the ith component of y at To pass the local error test, 
equation (2.91), the weighted local error vector, that is, {dij/EWT/i }, must 
satisfy the inequality 


1 N 

( 1 \ 2 

EWT., 


1 J 


(3.2) 


where EWT, j is the ith component of the error weight vector for the first step (see 
eq. (2.90)): ’ 


EWT- j = RTOL / |}< 0 | + ATOLj-. (3.3) 

In this equation RTOL/ and ATOL; are, respectively, the user- supplied local 
relative and absolute error tolerances for the ith solution component, F/o is the ith 
solution component at ^o, and the vertical bars |*| denote absolute value. 

The test given by equation (3.2) cannot be applied at the start of the step [%o, ] 

because y(£i) is not known. We therefore modify this test by using y(^o) as 
follows: We first define a weighted principal error function at order 1, <j>, with 
element given by 


where 


♦f- 


1 

2 W' ’ 


W i = EWT. ,/TOL, 


(3.4) 


(3.5) 
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Figure 3.2. — Flowchart of subroutine LSODE. 
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3. Description of Code 

and the scalar tolerance quantity TOL, which is to be determined, is such that Wj is 
a suitable weight for Y,, the ith component of X. The step size and the local error 
are then together required to satisfy the inequality 

A*|y<TOL, (3.6) 

where || # || represents a suitable norm. We have used a different symbol for the 
initial step size than in equation (3.2) to indicate that this quantity is not known 
and must be computed. Because a first-order method will be used on this step, for 
a sufficiently small step size the numerical approximation Xi at %\ will not be 
significantly different from y(J^o), and use of the latter quantity is therefore 
reasonable. The rationale for Introducing TOL will become apparent shortly. 

The second derivative y(^o) is not generally available, and so the following 
empirical procedure is used to estimate it. We consider the dominant eigenvalue 
(= X) of the ODE system and model this component with the simple scalar ODE 

y = = H ( 3 - 7 ) 

where | X | » 1. For this problem, <() = (1/2 )y/W = (U2)X 2 yfW. Now, if TOL is 
chosen such that y/W is of order unity, <)> can be approximated by (y/W) 2 
[= (Xy/W) 2 ], which is known. For the scalar ODE this condition is obtained by 
setting TOL = RTOL and ATOL = 0 (see eqs. (3.3) and (3.5)). The quantity y/W 
may be regarded as the weighted principal error function for a “zeroth order” 
method. We use this empirical rule to replace each $/ by ( y/W/) 2 so that equation 
(3.6) can be written as 


h l 



-|jV! 

-1 


< TOL, 


(3.8) 


where fj# [=fi(X 0, £o)l is the first derivative of the ith component at Because 
the weighted root-mean-square (rms) norm is used in the local error test, equa- 
tion (3.2), for convenience, we use the following criterion for initial step size 
control: 




i N 

ly 


ff V 2 

ho 

K W iJ 


< TOL. 


(3-9) 


Equations (3.5) and (3.9) together show that h 0 (<* 1/VTOL) is a decreasing 
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3.4 Special Features 

function of TOL. To produce a reliable estimate for h 0 , we therefore select a TOL 
erring on the high side. A suitable value is given by 


TOL = max(RTOL,). (3.10) 

i 


This expression cannot be used if all RTOL, = 0. In this case an appropriate value 
for TOL is given by 


TOL = max 


ATOL, 


V 


K 


for ^ 0 * 0. 


*,0 


(3.H) 


In any case the value of TOL is constrained to be within reasonable bounds as 
follows: 


100 u <. TOL < 1(T 3 , (3.12) 

where u is the unit roundoff of the computer or the machine epsilon (ref. 13). It is 
the smallest positive number such that 1 + n > 1. 

Equation (3.9) cannot be used to compute Hq if either each /, o is equal to zero 
or the norm is very small. To produce a reasonable ho in such an event, we 
include the independent variable C, as the zeroth component yo °f y an ^ modify 
equation (3.9) as follows: 


W r 


1 1 
b 

2 N 


-X 

hr t-t 


o 


N f f \ 

ho 

w. 

i = l V ' J 


< TOL, 


(3.13) 


where we have used the fact that yo = 1. To be consistent with the other Wj, 
which are of order y (j0 , the weight W 0 should be of order ^o; however, we use 

W 0 =max(|^o|^ou,,i|) < 314) 

to ensure that it is not equal to zero. In equation (3.14), ^out f l I s cither the first (or 
only) value of the independent variable at which the solution is required or, as 
discussed in chapter 4, a value that gives both the direction of integration (i.e., 
increasing or decreasing ^) and an approximate scale of the problem. If the 
quantity !^ u t i - is not significantly different from zero, an error exit occurs. 
Equation (3. i 3) gives a reasonable value for h 0 (= Wo VTOL) if fo = 0. 
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3. Description of Code 

The calculation procedure used for h§ is therefore given by 


h 0 = 


TOL 


%- 2 + 1 


vfA»l 


N' 


/ = ! 


W; 


' J 


(3.15) 


Several restrictions apply to the step size given by equation (3.15). It is not 
allowed to be greater than the difference |^o U t,l - £oi- Hence 

^<-min(v|5 ouM -M). (116) 

In addition, if the user has supplied a value for fc max , the maximum step size to be 
used on any step, Hq is restricted to 

fy ) *-min(V / w). ( 31? ) 

However, no comparison of Hq is made with h m j n , the user-supplied minimum 
step size to be used on any step, so that h 0 is allowed to be less than h m j n . Finally 
the sign of ho adjusted to reflect the direction of integration. 

3.4*2 Switching Methods 

Another useful feature of LSODE is that different integration methods and/or 
different iteration techniques can be used in different subintervals of the problem. 
This option is useful when the problem changes character and is stiff in some 
regimes and nonstiff in others as, for example, in combustion chemistry. Indeed, 
because stiff problems are usually characterized by a nonstiff initial “transient” 
region, the ability to switch integration methods is a desirable feature of any ODE 
package. During the course of solving a problem the method flag MF may be 
changed both whenever and as many times as desired. As described in chapter 4 
changing methods is quite straightforward. 

3.4.3 Excessive Accuracy Specification Test 

At each integration step LSODE checks that the user has not requested 

too much accuracy for the precision of the machine. This condition is said to 
occur if the criterion 

? 

d i n < uY. n (3.18) 

is true for all N solution components. In equation (3.18), di n is the estimated 
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3.4 Special Features 

local truncation error in Y, the ith solution component at ^ n . Now the numerical 
solution Xn at is judged to be sufficiently accurate if the following inequality is 
satisfied (see chapter 2): 



i N ( 

N &{ 


i.n 


EWT 


i.n ) 


? 

< 1 . 


(2.91) 


The quantity EWT/„ is the ith component of the error weight vector, equa- 
tion (2.90), for this step. Equations (3.18) and (2.91) together imply that if the 
quantity TOLSF (tolerance scale factor) defined as 


TOLSF = u. 


i N f 

-Y 

N ii{ 


i,n 


V EWT ,.ny 


(3.19) 


is greater than 1 , the test for excessive accuracy requirement is passed. This test is 
quite inexpensive, but it can be applied only after the solution at is produced. It 
is, however, wasteful to generate a solution only to discover that excessive 
accuracy has been required, either because TOLSF is greater than 1 or because 
repeated convergence failures or error test failures occur. The computational cost 
can be significant if any difficulty is encountered because of the corrective 
actions — described later in this section — performed by the code. Even if the step 
is successful, the solution is not meaningful because of roundoff errors. 

To avoid these difficulties, the calculation procedure for TOLSF uses Xi-b 
which is known, so that the test can be applied at the start of each step, including 
the first. Thus the code ascertains inexpensively if excessive accuracy has been 
requested before attempting to advance the solution by the next integration step. 
The value of TOLSF may be used to adjust the local error tolerances so that this 
condition does not recur. For example, scaling up the {RTOL,} and {ATOL,} 
values by a minimum factor of TOLSF should produce satisfactory values for the 
local error tolerances if the same type of error control is to be performed (see 
chapter 4 for details). 

3.4.4 Calculation of Method Coefficients 

The integration method coefficients and test constants used to check corrector 
convergence and local accuracy, as well as to select method order and step size, 
are computed in subroutine CFODE. The calculation procedure uses the generating 
polynomials discussed by Hindmarsh (refs. 21 and 22) to increase portability of 
the code. The coefficients corresponding to all method orders are computed and 
stored both at the start of the problem and whenever the user changes the 
integration method. This feature avoids the computational cost associated with 
recomputing these quantities whenever the method order is changed. 
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3. Description of Code 
3.4.5 Numerical Jacobians 


If Newton-Raphson (NR) or Jacobi-Newton (JN) iteration is selected, the code 
will generate elements of the Jacobian matrix by finite-difference approximations 
if the user chooses not to provide an analytical Jacobian. For the iteration 
procedures corresponding to MITER = 2 (full Jacobian matrix) and 5 (“banded” 
Jacobian matrix, i.e., a matrix with many zero entries and all nonzero elements 
concentrated near the main diagonal), the element Jy (= dfjfdyj) at is estimated 
by using the approximation 


Jus- 


fi + 5 kj AY j } .^] -Si ({ifi } ,5„] 


AY j 


* = 1,...,N, 


(3.20) 


where F^ rt is the kth component of Y^\ 8^ is the Kronecker symbol. 



k * j 
k = j* 


(3.21) 


and the increment A Yj in the jth solution component is selected as follows: The 
standard choice for A Yj is 


A Yj = 4~u\ 


y[0] 

j>n 


(3.22) 


This equation cannot be used if F^J is either equal to zero or very small. 
Therefore an alternative value, based on noise level, is deduced as follows: Now 
the error in each f due to roundoff is of order u\f\. Hence in replacing df/dyj by 
the difference quotient, equation (3.20), the resulting element /y has an error of 
order wj/jl/ry, where for clarity in presentation we have replaced AYj by rj. Finally 
because the method coefficient Jio (= <?o) is of order unity (see tables 2.1 and 2.2), 
the error 5 P x j in the element Py of the iteration matrix P, equation (2.25), is 
approximately 


5p u “ H“ \ f i\b (3 - 23) 

If we introduce the ^-dimensional column vector £, with element sj defined as 

Sj= \fr p ; = (3.24) 

the matrix 5P containing the errors {5 Py } is given by 
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5P = /z| u |f^ r , 


3.4 Special Features 


(3.25) 

where 1 f I is an Af-dimensional column vector containing the absolute values of the 
fi ( i = 1,...^/) and the superscript T indicates transpose. A suitable increment rj is 
obtained by bounding |8P|, as discussed next. 

To be consistent with the corrector convergence test, equation (2.98), and the 
local error test, equation (2.91), we use the weighted rms norm, which for an 
arbitrary /'/-dimensional column vector & is given by 


x = 




N ( \ 

X i 


EWT; 


(3.26) 


U 


If we introduce the diagonal matrix D of order N, with element Da given by 


D {i ~ 1/EWT-, i = A/, (3.27) 

it is easily verified that 

IWHNU^, (3.28) 

where ||*||£ is the Euclidean norm, defined for \ as 



(3.29) 


Now the norm of 5P is given by 


||SP|| = max 



where 


M - 


because 5P is of rank 1 . Hence 


(3.30) 


(3-31) 
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which can be rewritten as 



(3.32) 


To establish the maximum allowable error in P, we consider the linear system 
= h, which is the form of the equation to be solved at each Newton iteration, 
equation (2.24). To first order, the error 8* in * due to the error 8P in P is given by 
(e.g., ref. 13) 



(3.33) 


The norm | P _1 1 is not known but is expected to be of order unity because P I,, 
the identity matrix of order N t when h — » 0 and P h$ oJ when h — » °° (see 
eq. (2.25)) Therefore, a reasonable strategy is to bound 1 8P | alone by selecting a 
suitably small value for the relative error that can be tolerated in the Newton 
correction vector. By using a value of 0. 1 percent for this error, we obtain from 
equations (3.32) and (3.33) 


min 


EWT, 


i J 


10 3 |/i| uN ||f|| = r Q . 


(3.34) 


For additional safety r 0 is reset to 1 if it is equal to zero. Finally the increment A Yj 
in the jth variable used to estimate the {Jy} is given by 
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A Yj = max 



r 0 EWT ; 


],n 


(3.35) 


For a full Jacobian matrix the above procedure will require (N + 1) derivative 
evaluations and can therefore become much more expensive than the use of an 
analytical Jacobian, especially for large N. Now f(Xj? ] , Z, n ) is required by the 
corrector (see eq. (2.36)), irrespective of the iteration technique. Hence the use of 
MITER = 2 requires the evaluation of only N additional derivatives. 

In generating the finite-difference banded Jacobian matrix (MITER = 5) the 
code exploits the bandedness of the matrix for efficiency. The number of additional 
derivative evaluations required to form the Jacobian matrix is only ML + MU + 1 , 
where ML and MU are, respectively, the lower and upper half-bandwidths of the 
Jacobian matrix. 

If JN iteration with MITER = 3 is used, the ^diagonal elements/,; (i = I..../V) 
are estimated by using the approximation 



which requires only one additional derivative evaluation. The increment AT, is 
selected as follows: Now equation (2.17) shows that if functional iteration were 
used, the correction - Yj 01 that would be obtained on the first iteration is 
equal to the quantity PogOd° 5 )» where the vector function g is given by equa- 
tion (2. 16). The increment vector AX is taken to be 10 percent of this correction: 

A^-O.iPog,^ 01 ], i = (3.37) 


Hence the diagonal matrix approximation, equation (3.36), resembles a directional 
derivative of £ taken in the same direction as the correction vector above. Also, 
this approximation gives the correct Jacobian if it is a constant diagonal matrix. If 


the magnitude of A F/ is less than 0. 1 «Po EWT), that is, if 
set equal to zero. 




< wEWT f , Jn is 


3.4.6 Solution of Linear System of Equations 

If NR iteration is used for the problem, a linear system of the form Pa = h must 
be solved for the correction vector a at each iteration (see eq. (2.24)). The linear 
algebra necessary to solve this equation is performed by the LU method (e.g., 
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refs. 5 and 36), rather than by explicitly inverting the iteration matrix, which will 
require prohibitive amounts of computer time (ref. 13). In the LU method the 
iteration matrix is factored into the product of two triangular matrices L and U. 
Solving equation (2.24) then requires the fairly simple solution of two triangular 
linear systems in succession. 

LSODE also includes special procedures for the LU-decomposition of the 
iteration matrix and the solution of equation (2.24) when the matrix is known to 
be banded. Compared to a full matrix, it is significantly less expensive to form a 
banded matrix, perform its LU-decomposition, and solve the linear system of 
equations (refs, 5, 25, 26, and 36). An important advantage of LU-decomposing a 
banded matrix over inverting it is that, besides being faster, the triangular factors 
L and U lie within nearly the same bands as the original matrix, whereas the 
inverse is a full matrix (ref. 36). This feature makes the computation of the 
correction vector significantly faster with the LU method than by premultiplying 
the right-hand side of equation (2.24) with the inverse of the matrix. 

If MITER = 3 is used for the problem, the resulting iteration matrix is diagonal 
(see eq. (3.36)). Its inverse can therefore be obtained trivially and is used to 
compute the corrections. 

3.4.7 Jacobian Matrix Update 

The difficulty with Newton-Raphson iteration is the computational cost 
associated with forming the Jacobian matrix and the linear algebra required to 
solve for the correction vector at each iteration. However, as discussed in chap- 
ter 2, the iteration matrix need not be very accurate. This fact is exploited to 
reduce the computational work associated with linear algebra by not updating P at 
every iteration. For additional savings it is updated only when the iteration does 
not converge. Hence the iteration matrix is only accurate enough for the solution 
to converge, and the same matrix may be used over several steps. It is also 
updated if three or more error test failures occur on any step. Now P may be 
altered if the coefficient h% is changed (see eq. (2.25)) because a new step size 
and/or method order is selected. In order to minimize convergence failures 
caused by an inaccurate P, the code updates P and performs its LU-decomposition 
(or inversion if MITER = 3) if h$o has changed by more than 30 percent since the 
last update of P. In addition, for MITER = 3, because P" 1 can be generated 
inexpensively, it is first modified to account for any change in /i|3o since its last 
update, before the corrections are computed. The reevaluation and LU- 
decomposition or inversion are also done whenever the user changes any input 
parameter required by the code. Finally the same P is used for a maximum 
number of 20 steps, after which it is reevaluated and LU-decomposed or inverted. 

3.4.8 Corrector Iteration Convergence and Corrective Actions 

Irrespective of the solution method and the corrector iteration technique, the 
maximum number of corrector iterations attempted on any step is set equal to 3, 
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based on experience that a larger number increases the computational cost without 
a corresponding increase in the probability of successful convergence (refs. 19, 
21, 22, and 25). In addition to performing the convergence test, equation (2.99), 
at each iteration, STODE examines the value of the convergence rate c m , equa- 
tion (2.102). If c m is greater than 1, the iteration is clearly not converging. 
STODE exploits this fact by abandoning the iteration if c m is greater than 2 after 
the second iteration. 

If convergence is not obtained because either (1) equation (2.99) is not satisfied 
after three iterations or (2) c m > 2 after the second iteration, the following 
corrective actions are taken: For NR and JN iterations, if P is not current, it is 
updated at y = and LU-decomposed or inverted, and the step is retried with 
the same step size. However, if either P is current or functional iteration is used, a 
counter of convergence failures on the current step is increased by 1, the step size 
is reduced by a factor of 4, and the solution is attempted with the new step size. 
The same corrective actions are taken in the event of a singular iteration matrix. 

This procedure is repeated until either convergence is obtained or the integration 
is abandoned because either (1)10 convergence falures have occurred or (2) the 
step size has been reduced below a user-supplied minimum value /i m j n . In the 
event of an error exit the index of the component with largest magnitude in the 
weighted local error vector is returned to the subprogram calling LSODE. 

3.4.9 Local Truncation Error Test and Corrective Actions 

After successful convergence STODE performs the local truncation error test, 
equation (2.96). If the error test fails, the step size is reduced and/or the method 
order is reduced by 1 by using the procedures outlined in section 3.4.10, and the 
step is retried. After two consecutive failures the step size is reduced by at least a 
factor of 5, and the step is retried with either the same or a reduced order. After 
three or more failures it is assumed that the derivatives that have accumulated in 
the Nordsieck history matrix have errors of the wrong order. Therefore the First 
derivative is recomputed and the method order is set equal to 1 if it is greater than 
1. Then the step size is reduced by a factor of 10, the iteration matrix is formed 
and either LU-decomposed or inverted, and the step is retried with a new z n „\ that 
is constructed from Xz-1 and Y fl -i = £QC*-i). 

This procedure is repeated until either the error test is passed or an error exit is 
taken because either (1)10 error test failures have occurred or (2) the step size has 
been reduced below In the event of an error exit LSODE returns the index of 
the component with the largest magnitude in the weighted local error vector to the 
calling subprogram. 

If the accuracy test is passed, the step is accepted as successful, and the 
Nordsieck history matrix z n and the estimated local truncation error vector d* at ^ 
are computed by using equations (2.76) and (2.89), respectively. Irrespective of 
whether the step was successful or not, STODE saves the value of the most recent 
step size attempted on the step so that the user may, if desired, change it. 
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3.4.10 Step Size and Method Order Selection 

In addition to advancing the solution STODE periodically computes the method 
order and step size that together maximize efficiency while maintaining prescribed 
accuracy. As discussed in chapter 2, this result is accomplished by selecting the 
method order that maximizes step size. To simplify the algorithm, the code 
considers only the three method orders q - 1, q, and q + 1, where q is the current 
method order. For each method order the step size that will satisfy exactly the 
local error bound is computed by assuming that the highest derivative remains 
constant. The resulting step size ratios (defined as the ratio of the step size to be 
attempted on the next step to the current value h n ) are given by equations (2. 107), 
(2.103), and (2.1 12), respectively, for method orders q - 1, q, and q + 1. These 
equations are, however, modified by using certain safety factors (1) to produce a 
smaller step size than the value that satisfies the error bound exactly, because the 
error estimates are not exact and the highest derivative is not usually constant, and 
(2) to bias the order-changing decision in favor of not changing the order at all, 
because any change in order requires additional work, and then in favor of 
decreasing the order, because an order reduction results in less work per subsequent 
step than an order increase. The formulas used in STODE to calculate the step 
size ratios are as follows: 


'down 


r up = 


1.3 

1 

( D q -i) q +10" 6 


1 

1.2 

1 

i +io-* 

r~ 

i 


1.4 


H + i) 9+2+10 ^ 


(3.38) 


(3.39) 


(3.40) 


In equations (3.38) to (3.40) the factors 1.2, 1.3, 1.4, and 10 -6 are strictly 
empirical. The subscripts “down,” “same,” and “up” indicate, respectively, that 
the method order is to be reduced by 1, left unchanged, and increased by 1. 

To prevent an order increase either after a failed step or when q = q max , the 
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maximum order allowed for the solution method, r up is set equal to zero in such 
cases. Similarly, if q = 1, r^own is set equal to zero to avoid an order reduction. 

The maximum step size ratio r = max (rdowm Tsame* r up) and the corresponding 
method order are selected to be attempted on the next step if r ^ 1.1 after a 
successful step. Changes in both step size and method order are rejected if the 
step size increase is less than 10 percent because it is not considered large enough 
to justify the computational cost required by either change (refs. 10 and 22). After 
a failed step the method order is decreased if r<] own > r same ; however, r - max 
( r down> r same) is reset to 1 if it is greater than 1. Several additional tests, given 
next, are performed on r, if r ^ 1 . 1 after a successful step, but irrespective of the 
value of r after a failed step, before the step size h' (= rh n ) to be attempted next is 
selected. 

If the maximum step size h max to be attempted on any step has been specified 
by the user, r is restricted to 


mini 


r, - 


(3.41) 


Similarly if the user has specified a minimum step size h m \ n that may be attempted 
on any step, r is restricted to 


r <— max 


r>- 


n n ) 


Finally r must satisfy the inequality 


r < r 


max, 


(3.42) 


(3.43) 


where the variable r max is normally set equal to 10. However, for the very first 
step size increase for the problem, if no convergence or error test failure has 
occurred, r max is set equal to 10 4 to compensate for the small step size attempted 
on the first step. For the first step size increase following either a corrector 
convergence failure or a truncation error test failure, r max is set equal to 2 to 
inhibit a recurrence of the failure. 

To avoid numerical instability caused by frequent changes in the step size, 
method order and step size changes are attempted only after S successful steps 
with the same method order and step size, where 5 is normally set equal to q + 1 . 
However, if an unsuccessful step occurs, this rule is disregarded and the step size 
and/or the method order may be reduced. Following a failed error test or a failed 
convergence test with either functional iteration or NR and JN iterations if P is 
current, S is set equal to q + 1 . If three or more error test failures occur on any one 
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step, S is set equal to 5 even though the method order is reduced to 1. Finally 
following a step for which step size and method order changes are rejected 
because r < 1 . 1 , S is set equal to 3. 

After every S - I successful steps STODE saves the vector e, if q < q max , in 
order to estimate Vs, which is required to compute r up (see eqs. (2.109) to 
(2.112)). To minimize storage requirements, £„ is saved as the 9 max th, that is, the 
last, column of z„. 


3.5 Error Messages 

The code contains many error messages — too numerous to list here. Every 
input parameter is tested for legality and consistency with the other input variables. 
If an illegal input parameter is discovered, a detailed message is printed. Each 
error message is self-explanatory and complete. It not only describes the mistake 
but in some instances tells the user how to fix the problem. Any difficulty 
encountered during execution will result in an error exit. A message giving the 
reason for termination will also be printed. If the computation stops prematurely, 
the user should look for the error message near the end of the output file 
corresponding to the logical unit number LUNIT (see chapter 4). 
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Description of Code Usage 

To use the LSODE package, the following subprograms must be provided: (1) a 
routine that manages the calls to subroutine LSODE, (2) a routine that computes 
the derivatives {// = dyjdt) for given values of the independent variable ^ and the 
solution vector y , and (3) if an analytical Jacobian matrix J (= di/d y) is required 
by the corrector iteration technique selected by the user, a routine that computes 
the elements of this matrix. In addition, some modifications, discussed below, to 
the LSODE source itself may be necessary. 


4.1 Code Installation 

4.1.1 BLOCK DATA Variables 

The user may wish to reset the values for the integer variables MESFLG (cur- 
rently 1 ) and LUNTT (currently 6), which are both set either in the BLOCK DATA 
module (double-precision version) or in subroutine XERRWV (single-precision 
version). The variable MESFLG controls the printing of error messages from the 
code, and LUNIT is the logical unit number for such output (see table 3,7). 
Setting MESFTG = 0 will switch off all output from the code and therefore is not 
recommended. 

The single-precision version of the code loads initial values for the common 
block LS0001 variables ILLIN and NTREP (see table 3.8) through a DATA state- 
ment in subroutine LSODE. The same procedure is used in subroutine XERRWV 
for the common block EH0001 variables MESFLG and LUNIT (see table 3.7). 
However, on some computer systems initial values for common block elements 
cannot be defined by means of DATA statements outside a BLOCK DATA 
subprogram. In this case the user must provide a separate BLOCK DATA 
subprogram, to which the two DATA statements from subroutines LSODE and 
XERRWV must be moved. The BLOCK DATA subprogram must also contain 
the two common blocks EH0001 and LS0001 (see table 3.6). 
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4.1.2 Modifying Subroutine XERRWY 

The subroutine XERRWV, which prints error messages from the code, is 
machine and language dependent. Therefore the data type declaration for the 
argument MSG, which is a Hollerith literal or integer array containing the message 
to be printed, may have to be changed. The number of Hollerith characters stored 
per word is assumed to be 4, and the value of NMES, which is the length of, that 
is, number of characters in, MSG is assumed to be a multiple of 4, and at most 60. 
However, the routine describes the necessary modifications for several machine 
environments. In particular, the user must change a DATA statement and the 
format of statement number 10. The routine assumes that all errors are either (1) 
recoverable, in which case control returns to the calling subprogram, or (2) fatal, 
in which case the run is aborted by passing control to the statement STOP, which 
may be machine dependent. If a different run-abort command is needed, the line 
following statement number 100, which is located near the end of the routine, 
must be changed. 


4.2 Call Sequence 

The call sequence to subroutine LSODE is as follows: 

CALL LSODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, ISTATE, 
IOPT, RWORK, LRW, IWORK, LIW, JAC, MF) 

All arguments in the call sequence are used on input, but only Y, T, ISTATE, 
RWORK, and IWORK are used on output. Also, Y and T are set only on the first 
call to LSODE; the other arguments may, however, have to be reset on subsequent 
calls. The arguments to LSODE are defined as follows: 

F The name of the user-supplied subroutine that computes the derivatives 

of the dependent variables with respect to the independent variable. 
This name must be declared EXTERNAL in the subprogram calling 
LSODE. The requirements of subroutine F are described in section 
4.3. 

NEQ The number of first-order ordinary differential equations (ODE’s) to 
be solved. (The code allows the user to decrease the value of NEQ 
during the course of solving the problem. This option is useful if 
some variables can be discarded as the solution evolves as, for example, 
in chemical kinetics problems for which the reaction mechanism is 
reduced dynamically.) As discussed later, NEQ can be specified as an 
array. In this case NEQ(l) must give the number of ODE’s to be 
solved, and the subprogram calling LSODE must contain a dimension 
statement for NEQ. 
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Y A vector of length NEQ (or more) containing the dependent variables. 

The subprogram calling LSODE must include a dimension statement 
for Y if it contains more than one component. On the first call to 
LSODE this vector must be set equal to the vector of initial values of 
the dependent variables. Upon every return from LSODE, Y is the 
solution vector either at the desired value (TOUT or TCRTT, see 
below) of the independent variable or that generated at the end of the 
previous integration step. In case of an error exit Y contains the 
solution at the last step successfully completed by the integrator. 

T The independent variable. On the first call to LSODE, T must give 

the initial value of this variable. On every return from LSODE, T is 
either the independent variable value (TOUT or TCRIT, see below) at 
which the solution is desired or the independent variable value to 
which the numerical solution was advanced on the previous integration 
step. If an error exit occurs, T gives the value of the farthest point (in 
the direction of integration) reached by the integrator. 

TOUT The next value of the independent variable at which the solution is 
required, if ITASK = 1, 3, or 4 (see table 4.1). For ITASK = 2 or 5, 
LSODE uses TOUT on the first call to determine the direction of 
integration and, if necessary, to compute the step size to be attempted 
on the first step; on subsequent calls TOUT is ignored. LSODE 
permits integration in either direction of the independent variable. 

ITOL A flag that indicates the type of local error control to be performed. 

The legal values that can be assigned for ITOL and their meanings are 


TABLE 4.1. — VALUES OF ITASK USED IN LSODE 
AND THEIR MEANINGS 


ITASK 

Description 

a l 

Compute output values of Y(^) at ^ by overshooting and 

interpolation. 

2 

Advance the solution to the ODE’s by one step and return to 
calling subprogram. 

a 3 

Stop at the first internal mesh point at or beyond £ - and 

return to calling subprogram. 

aJ>4 

Compute output values of Y<£) at £ = but without over- 

shooting \ = ^ ri( . 

b 5 

Advance the solution to the ODE's by one step without passing 
% - and return to calling subprogram. 


a User must supply value for (= TOUT). 

kUser must supply value for ^ rit (= TCRIT)- This option is useful if the 
problem has a singularity at or beyond £ = ^ rit . 
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RTOL 


ATOL 


ITASK 


ISTATE 


TABLE 4.2.— VALUES OF ITOL USED 
IN LSODE AND THEIR MEANINGS 


ITOL 

Description 

1 

Scalar RTOL and scalar ATOL 

2 

Scalar RTOL and array ATOL 

3 

Array RTOL and scalar ATOL 

4 

Array RTOL and array ATOL 


given in table 4.2. The variables RTOL and ATOL are described next. 

The local relative error tolerance parameter for the solution. This param- 
eter can be specified either as a scalar, so that the same tolerance is used 
for all dependent variables, or as any array of length NEQ, so that 
different tolerances are used for different variables. In the latter case the 
subprogram calling LSODE must contain a dimension statement for 
RTOL. 

The local absolute error tolerance parameter for the solution. This 
parameter can also be specified either as a scalar, so that the same 
tolerance is used for all dependent variables, or as an array of length 
NEQ, so that different tolerances are used for different variables. In 
the latter case the subprogram calling LSODE must contain a dimension 
statement for ATOL. 

An index that specifies the task to be performed. This flag controls 
when LSODE stops the integration and returns the solution to the 
calling subprogram. The legal values for ITASK and their meanings 
are given in table 4.1. If ITASK = 4 or 5, the input variable TCRIT (= 
independent variable value that the integrator must not overshoot, see 
table 4.1) must be passed to LSODE as the first element of the array 
RWORK (defined below). 

An index that specifies the state of the calculation, that is, if the call to 
LSODE is the first one for the problem or if it is a continuation. The 
legal values for ISTATE that can be used on input and their meanings 
are given in table 4.3. The option ISTATE = 3 allows changes in the 
input parameters NEQ, ITOL, RTOL, ATOL, IOPT, MF, ML, and MU 
and any optional input parameter, except HO, discussed in the 
descriptions of RWORK and IWORK. The integer variables IOPT, 
MF, ML, and MU are defined below. The parameters ITOL, RTOL, 
and ATOL may also be changed with ISTATE = 2, but LSODE does 
not then check the legality of the new values. On return from LSODE, 
ISTATE has the values and meanings given in table 4.4. 
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ISTATE 

Description 

1 

2 

3 

ar , 

This is the first call for the problem 

TOs is not the first call for the problem, and the calculation is to 
be continued normally with no change in any input paralters 
1 exce P< Possibly ^ and ITASK. 4 parameters 

This is not the first call for the problem, and the calculation is to 
be continual normally, but with a change in input parameters 

other than ^ and ITASK 4 


See table 4.1 for description of IT ASK. 


TABLE 4.4 — VALUES OF ISTATE RETURNED BY LSODE 
and THEIR MEANINGS 


ISTATE 


2 

-1 


-2 

-3 

-4 

“5 

~6 


Meaning 


a See table 4.6. 
b See chapter 2. 


Noting was done because TOUT = T on first call to LSODE 
(However, an internal counter was set to detect and prevent 
repeated calls of this type.) 

The integration was performed successfully. 

Excessive amount of work was done on this call (i.e., number of 
eps exceeded MXSTEP 4 on this call), bu, the integration w* 
successful as far as the value returned in T 

Too much accuracy was requested for the computer being used but 
the integrate was successful as far as the value returned in T. 

IS error is detected on the first call to LSODE (i e before 

*' " i "' g ' 1 *" - 3 - ~ 

Repeated eiror test failures occurred on one step, but the integration 
was successful as far as the value returned in V g 

Repeated convergence test failures occurred on one step, but the 
integration was successful as far as the value returned in T 
Some component, EWT„ of the error weight vector EWT* ' 
vanished, so that the local error test cannot be app^Tbut the 
integration was successful as far as the value retied in T (This 
condmon arises when pure relative error control (i.e.. ATOL 
- -2 WaS Specifled fof a vanabl <= whose magnitude is now zero ) 
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:r r '^rrrL g that specifics if any optional input is being used on 
to S The legal values for IOPT together with the,, meanmgstue 
„jven in table 4.5. The optional input parameters that may be set by 
*e use, are given tn table 4.6. For each such input variable this table 
lists its location in the call sequence, its meaning, and its 
late. The quantities RWORK and IWORK arc work anays descnbed 

below. 


TABLE 4 5 -VALUES OF IOPT THAT CAN RE USED ON 
TABLE, tuctd XfFANINGS 


INPUT TO LSUUfc 1 

IOPT 

Meaning 

0 

1 

The user has not set a value for any optional input parameter* 
(Default values will be used for all these parameters.) 
VtoTL. tea speciM ,» - « 

parameters. - 


a See table 4.6 for a list of these parameters. 


TABLE 4.6. 
AND 


AtmtlNAl INPUT PARAMETERS THAT CAN BE SET BY USER 
ijSSoxEg tSS «« default VALUES, 


Optional 

input 

parameter 


HO 

HMAX 

HMIN 

MAXORD 

MX STEP 
MXHNIL 


Location 


RWORK(5) 

RWORK(6) 

RWORK(7) 

IW0RK(5) 

lWORK(6) 

IW0RK(7) 


Meaning 


Step size to be attempted on 
the first step 

Absolute value of largest step 
size (in magnitude) to be 
used on any step 
Absolute value of smallest 
step size (in magnitude) to 
be used on any step a 
Maximum method order to be 
used on any step 

Maximum number of integra- 
tion steps allowed on any 
one call to LSODE 
Maximum number of times 
that warning message that 
step size is getting too small 
is printed 


Default value 


Computed by LSODE 


12 for Adams-Moulton 
method and 5 for 
backward differenti- 
ation formula method 
500 


10 


•m. ,*» U ignted oa the M «P te <* K “ •* “ ™ Wl “ 

XXASK = 4 or 5 (see table 4.1). 
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RWORK A real work array used by the integrator. The subprogram calling 
LSODE must include a dimension statement for RWORK. If ITASK = 
4 or 5, the user must set RWORK(l) = TCRIT (see table 4.1) to 
transmit this variable to LSODE. If any optional real input parameters 
are used, their values are also passed in this array to LSODE; the 
address for each of these parameters is given in table 4.6. Upon return 
from LSODE, RWORK contains several optional real output 
parameters. For each such output variable table 4.7 lists its location in 
RWORK and its meaning. In addition, the Nordsieck history array at 
the current value of the independent variable (TCUR in table 4.7) and 
the estimated local error vector in the solution incurred on the last 
successful step can be obtained from RWORK. Table 4.8 lists the 
names used for these two quantities and their locations in RWORK. 
In this table NYH is the value of NEQ on the first call to LSODE, and 
NQCUR and LENRW are both defined in table 4.7, which also gives 
their locations in the array IWORK (see below). 

LRW Length of the real work array RWORK. Its minimum value depends 
on the method flag MF (see below) and is given in table 4.9 for each 
legal value of MF. In this table the integer MAXORD is the maximum 
method order (default values = 1 2 and 5 for the AM and BDF methods, 
respectively) to be used. The integers ML and MU are the lower and 
upper half-bandwidths, respectively, of the Jacobian matrix if it is 
declared to be banded (see table 3.2). 

IWORK An integer work array used by the integrator. The subprogram calling 
LSODE must include a dimension statement for IWORK. If MITER 
(— second decimal digit of MF, defined below) = 4 or 5 (table 3 2) the 
user must set IWORK( 1) = ML and IWORK(2) = MU (see descriptions 
above) to transmit these variables to LSODE. If any optional integer 
input parameters are used, their values are also passed in this array to 
LSODE; the address for each of these parameters is given in table 4.6. 
Upon return from LSODE, IWORK contains several optional integer 
output parameters. For each such output variable table 4.7 lists its 
location in IWORK and its meaning. 

LIW Length of the integer work array IWORK. Its minimum value depends 

on MITER (table 3.2) and is given in table 4. 10 for each legal value of 
MITER. 

AC The name of the user-supplied subroutine that computes the elements 
of the Jacobian matrix. This name must be declared EXTERNAL in 
the subprogram calling LSODE. The form and description of sub- 
routine JAC are given in section 4.4, 
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TABLE 4.7.— OPTIONAL OUTPUT PARAMETERS RETURNED BY LSODE 
AND THEIR LOCATIONS AND MEANINGS 


Optional 

output 

parameter 

Location 

Meaning 

HU 

RWORK(ll) 

Step size used on last successful step 

HCUR 

RWORK(12) 

Step size to be attempted on next step 

TCUR 

RWORK(13) 

Current value of independent variable. The 
integrator has successfully advanced the 
solution to this point. 

TOLSF 

RWORK(14) 

A tolerance scale factor, greater than 1 .0, that 
is computed when too much accuracy is 
requested (1ST ATE = -2 or -3, see table 4.4). 
To continue integration with the same ITOL, 
the local error tolerance parameters RTOL 
and ATOL must both be increased by at 
least a factor of TOLSF. 

NST 

IWORK(ll) ] 

Number of integration steps used so far for 

NFE 

IWORK(12) 

problem 

Number of derivative evaluations required so 
far for problem 

NJE 

IWORK(l3) 

Number of Jacobian matrix evaluations (and 
iteration matrix LU -decompositions or 
inversions) so far for problem 

NQU 

IWORK(14) 

Method order used on last successful step 

NQCUR 

IWQRK(15) 

Method order to be attempted on next step 

IMXER 

IWORK(16) 

Index of component with largest magnitude in 
weighted local error vector (e/EWT,, see 
chapter 2). This quantity is computed when 
repeated convergence or local error test 
failures occur. 

LENRW 

IWORK(17) 

Required length for array RWORK 

LENTW 

IWORK(18) 

Required length for array IWORK 


TABLE 4,8. USEFUL INFORMATIONAL QUANTITIES REGARDING INTEGRATION 

THAT CAN BE OBTAINED FROM ARRAY RWORK 


Quantity 

Name 

Location 

Nordsieck history array for problem 

Estimated local error in solution on 
last successful step 

YH 

ACOR 

RWORK(2!) to 

RWORK(20 + NYH(NQCUR + 1» 
RWORK(LENRW - NEQ + 1) to 
RWORK (LENRW) 
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4.3 User-Supplied Subroutine for Derivatives (F) 
TABLE 4.9 — MINIMUM LENGTH REQUIRED BY REAL WORK 


MF 

Minimum LRW® 

10,20 

11,12,21,22 

13,23 

14,15,24,25 

20 + NYH(MAXORD + 1) + 3 NEQ 

22 + NYH(MAXORD + 1) + 3 NEQ + (NEQ) 2 

22 + NYH(MAXORD + 1) + 4 NEQ 

22 + NYH(MAXORD + 1) + (2 ML + MU + 4)NEQ 


uh urei can to LoUDE. 

MAXORD is the maximum method order to be used for problem 
NEQ is the number of ODE's specified on current call to LSODE, 
and ML and MU are. respectively, the lower and upper haif- 
b and widths of the banded Jacobian matrix. 


MF 


TABLE 4.10.— MINIMUM 
LENGTH REQUIRED BY 
INTEGER WORK ARRAY 
IWORK (i.e., MINIMUM 
LIW) FOR EACH MITER 


MITER* 

Minimum LIW b 

0 

20 

1,2 

20 + NEQ 

3 

20 

4,5 

20 + NEQ 


“See table 3.2 for description 
of MITER. 

b NEQ is the number of ODE’s 
specified on current call to 
LSODE. 


Method flag that indicates both the integration method and corrector 
Ration technique to be used. MF consists of the two decimal digits 
MbT ”’ whlch s P ec 'fies the integration method, and MITER, which 
specifies the iteration technique (eq. (3.1)). Equation (3.1) and 

in’ n n Tf II . S J°^ that MF has the Allowing 12 legal values— 

, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, and 25. If MF= 14, 15 24 or 

25 the values of ML and MU must be passed to LSODE as the first 
and second elements, respectively, of the array IWORK (see above) 


4.3 User-Supplied Subroutine for Derivatives (F) 


solve the nmN Sf S ° °" method ° r corrector itera(i on technique selected to 
olve the problem the user must provide a subroutine that computes the derivatives 

^ T ° f thC independent var iaMe and the solution vector. The 
(F) o this subroutine is an argument in the call vector to LSODE and must 
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EXT^AL in .he subprogram calling LSODE. The 
derivative subroutine F should have the form 

SUBROUTINE F (NEQ, T, Y, YD OJ) 

statement for it. The routi where N is the current number of 

with YDOT(I) -dyjdt, (i I). evaluatedat ’j^ nt j t j es w hose current values, 

If the calculation of externa lly to LSODE, a special calculation, 

that is, at (or South *“* 9 , X h e results of the last call from the 

such as a cal! to .he rouhne . F, mu* * ' corrKp o„d a Y value 
package to the routine F shoe b a , va]ue (ha| may different from 

rd“i^ iSMen. variable value .o whmh .h^nmerica. 

.solution was advanced on .Ins previous ^o^ok,™ ihTsMragt requirement, the 
special call to subroutine F is ma , rwoRK(LSAVF) the base address of an 
YDOT argument ^ IJSSSSSU »> 

N-dimensional array, SAvF ( words) in the common 

LSAVF is the 224th word (6.h .meger word 
block LS0001 (table 3.6). If the derivative Y„ is required, 
calling subroutine INTDY, as explained in section 4.8. 

4.4 User-Supplied Subroutine for Analytical 
Jacobian (JAC) 

If .he corrector deration technique seized by ^ ^ ^ajlwtobian 

matrix, we ^“etlAC) o/thisroutine is an argument in the call vector to 
LSODlf and must'dmrefore be declared EXTERNAL in the snbprogmm call.ng 
LSODE. The Jacobian subroutine JAC should have t e orm 

SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD^^ D > 
DIMENSION Y(l), PD (NROWPD, 1) in FORTIN 66 
DIMENSION Y(»), PD (NROWPD, *) m FORTRAN 

Here ML and MU are, respectively; the by 
bandwidths of the Jacobian matnx if it is banded, and NKUW 
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Srottod f MU F " n b a„ n d e f d m Til 

™'V f NEQ is “ 

altered. For a full Jacobian mania (MITER = 1) the element rocW^’T' £ 

J = 1.....N) must be loaded with fyfiyL T v ( t = I .,_j\ T ... ’" ’ 

1 -\/|S=T;y=Y\‘ In this case the 

arguments ML and MU are not nepd^rl Tf t u* 

- 4^ th* i -w ^ d d ‘ If he Jacobian matrix is banded (MITER 

loaded in column-wise manned with^jTlt « tf ' ’j "fr” T ™ T " * 
loaded into the tows of PD. For a diagonafmani ' l U - old ^ 

elements must be loaded into a single row of length w t ° °’ dia S onal 

all elements of PD equal to zero before calling MC soThat^Iv ** 

elements need to be loaded Also eaoh ™i 1 1 if ■ th 1 y the nonzer ° 
.O subroutine F with “ Pre “*‘' by 8 

efficiency, intermediate quantises needed 

method. accessed by JAC by means of this 

^nchon.' iteration ^R = 0) or an internally generated Jacobian matrix 

RETURN TINE JAC (NEQ ’ T ’ Y ’ ML ’ MU ’ P °’ NROW PD) 

END 


4.5 Detailed Usage Notes 

It is apparent from the description of the call sequence to IS OOF ,w m j 

~;S~S=SS 

(ref. 37), who studied the effects ir P b em ‘yi* 5 ’ and by Radhalcrishnan 
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4. Description of Code Usage 
4.5.1 Normal Usage Mode 

The normal mode of communion with LSODE may be summaneed as 
follows: 

(1) Set initial values in Y. > and MF. 

(2) Set NEQ, T, UOL. RTOL, AT°ULRW,U = 1>and IOPT = 0 . 

(3) Set TOUT = first output station, 1TASK 

(4) Call LSODE. 

(5) Exit if 1STATE < 0. 

(6) Do desired output of Y. 

$ SSSSS “-on and return to srep (4). 

This procedure will result in LSODE (a) 

on the first step, (b) continuing the Jond TOUT, and (c) 
until the first internal meshponn ■ The returned value T will be 

computing the solution at TOUT by solution a t TOUT. Because the 

set equal to TOUT exac ly andY w.llc ^ ^ ^ ^ ^ for normal 

normal output value of IS! Alt, is a 
continuation. 

4.5.2 Use of Other Options 

1 m aif P use of other options included in the 
The calling subprogram may ' a Is ^ ^ ^ reS et to 3 to indicate that at 

package. For example, m step have ^ chang ed. The task to be 

TOUT some parameters, such as N Q - ’ Qm however , be changed without 

performed, indicated by the va u ® ' “ d ’ difficulties parameter values may 

prevent a tecurmnce of Ore indicated trouble. 


4.5.3 Dimensioning Variables 


UlllIV*»ur. 

Irrespective of the options selected^ the ^encT variables 8 that are arrays, 
include DIMENSION statemen s RWORK, IWORK, and, as discussed 

Such variables include Y, RT > > declared to be of length NEQ 

below, possibly NEQ. The ^n^orY^ de whose 

— e 

?££££ dements of Y; the remaining eiements are unchanged 

^ The parameter NEQ is nsuaii, a scaiar puantity. However, an array NEQ may 
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4.5 Detailed Usage Notes 

be used to store and pass integer data to the routines F and/or JAC. In this case the 
first element of NEQ must be set equal to the number of ODE’s. The LSODE 
package accesses only NEQ(l). However, NEQ is used as an argument in the 
calls to the routines F and JAC, so that these routines, and the MAIN program, 
must include NEQ in a DIMENSION statement. 

4.5.4 Decreasing the Number of Differential Equations (NEQ) 

In the course of solving a problem the user may decrease (but not increase) the 
number of ODE’s. This option is useful if some variables reach steady-state 
values while others are still varying. Dropping these constant quantities from the 
ODE list decreases the size of the system and hence increases computational 
efficiency. To use this option, upon return from LSODE at the appropriate time, 
the calling subprogram must reset the value of NEQ (orNEQ(l)); set ISTATE = 3; 
reset the values of all other parameters that are either required to continue the 
integration, such as TOUT if ITASK = 1 , 3, or 4 (table 4. 1 ), or are changed at the 
user’s option; and then call LSODE again. If the Jacobian matrixes declared to be 
banded (MITER = 4 or 5, table 3.2) and reductions can be made to the half- 
bandwidths ML and MU, they will also produce efficiency increases. The option 
of decreasing the number of ODE’s may be exercised as often as the user wishes. 
Of course, each time the size of the ODE system is decreased the changes 
discussed above should be made and the resulting number of ODE’s can never be 
less than 1 . However, the LRW and LIW values need not be reset. 

If, at any time, the number of ODE’s is decreased from N to N', LSODE will 
drop the last N- N'ODE’s from the system and integrate the first N’equations. It 
is therefore important in formulating the problem to order the variables carefully 
and make sure that it is indeed the last N - N' variables that attain steady-state 
values. In continuing the integration LSODE will access only the first N' elements 
of Y. However, the remaining N — N', or more, elements can be accessed by the 
user, and so no special programming is needed in either routine For JAC. 

4.5.5 Specification of Output Station (TOUT) 

The argument TOUT must be reset every time LSODE is called if the option 
given by ITASK = 1, 3, or 4 is selected. For the other two values of ITASK (i.e., 2 
and 5), TOUT need be set only on the first call to LSODE. Irrespective of the 
value of ITASK, the TOUT value provided on the first call to LSODE is used to 
determine the direction of integration and, if the user has not supplied a value for 
it, to compute the step size to be attempted on the first step. Therefore unless the 
user specifies the value for the initial step size, it is recommended that some 
thought be given to the value used for TOUT on the first call to LSODE. 

On the first call to LSODE, that is, with ISTATE = 1, TOUT may be set equal to 
the initial value of the independent variable. In this case LSODE will do nothing, 
and so the value ISTATE = 1 will be returned to the calling subprogram; however, 
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4. Description of Code Usage 

an internal counter will be updated to prevent repeated calls of this nature. If such 
a “first” call is made more than four times in a row, an error message will be 

issued and the execution terminated. . , 

On the second and subsequent calls to LSODE there is no requirement that the 
TOUT values be monotonic. However, a value for TOUT t at ac si up is 
limited to the current internal interval [(TCUR - HU),TCUR], where TOUR is he 
current value of the independent variable and HU is the step size used on the 

previous step. 

4,5.6 Specification of Critical Stopping Point (TOUT) 

In addition to TOUT a value must be specified for TCRIT if the option 
ITASK = 4 is selected. TCRIT may be equal to TOUT or beyond it, but not 
behind it, in the direction of integration. The integration is not permitted to 
overshoot TCRIT, so that the option is useful if, for example, a sin S u '^ y exi 
at or beyond TCRIT. This variable is also required with the option ITASK. — j. in 
either case the first element of the array RWORK (i.e., RWORK(l)) must be : set 
equal to TCRIT. If the solver reaches TCRIT within roundoff it will retu 
T = TCRIT exactly and the solution at TCRIT is returned in Y. o continue 
integrating beyond TCRIT, the user must reset either ITASK or TCRIT. In either 
case the value of ISTATE need not be reset. However, whenever TCRIT 
changed, the new value must be loaded into RWORK(l). 

4.5.7 Selection of Local Error Control Parameters (ITOL, RTOL, and 
ATOL) 

Careful thought should be given to the choice of ITOL, which together ^wjth 
RTOL and ATOL determines the nature of the error control performed by L • 

The value of ITOL dictates the value of the local error weight vector E fflCL with 
element EWT,- defined as 


EWT. = RTOL j |^| + ATOL ( , 


(4.1) 


where RTOL; and ATOL; are, respectively, the local relative and absolute e 
tolerances for the ith solution component Y, and the bars |-| denote absolute value. 
The solver controls the estimated local errors {d,} in {T;} by requiring the root- 

mean-square (rms) norm of <f, /EWT, to be 1 or less. . 

Le relative error control for the ith solution component ts obtained b, setting 
ATOL = 0- RTOL, is then a measure of the number of accurate sigmfican lg 
ures in The ’numerical solution. This error control is generally appropriate when 

widely varying orders of magnitude in 7, are expected. However, it cannot be 

used if the solution vanishes because relative error is then undefine . 
absolute error control for the ith solution component is obtained by setting 
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4.5 Detailed Usage Notes 

RTOL, = 0; ATOL, is then a measure of the largest number that may be neglected. 

Both RTOL and ATOL Can be specified (1) as scalars, so that the same error 
tolerances are used for all variables, or (2) as arrays, so that different tolerances 
are used for different variables. The value of the user-supplied parameter ITOL 
indicates whether RTOL and ATOL are scalars or arrays. The legal values that 
can be assigned to ITOL and the corresponding types of RTOL and ATOL are 
given in table 4.2. If RTOL and/or ATOL are arrays, the calling subprogram must 
include an appropriate DIMENSION statement. A scalar RTOL is generally 
appropriate if the same number of significant figures is acceptable for all 
components of Y. A scalar ATOL is generally appropriate when all components of 
Y, or at least their peak values, are expected to be of the same magnitude. 

In addition to ITOL, RTOL and ATOL should be selected with care. Now the 
code controls an estimate of only the local error, that is, an estimate of the error 
committed on taking a single step, starting with data regarded as exact. However, 
what is of interest to the user is the global truncation error or the actual deviation 
of the numerical solution from the exact solution. This error accumulates in a 
nontrivial manner from the local errors and is neither measured nor controlled by 
the code. It is therefore recommended that the user be conservative in choosing 
values for the local error tolerance parameters. However, requesting too much 
accuracy for the precision of the machine will result in an error exit (table 4.4). In 
such an event the minimum factor TOLSF by which RTOL and ATOL should both 
be scaled up is returned by LSODE (see table 4.7). Some experimentation may be 
necessary to optimize the tolerance parameters, that is, to determine values that 
produce sufficiently accurate solutions while minimizing the execution time. The 
global errors in solutions generated with particular values for the local error 
tolerance parameters can be estimated by comparing them with results produced 
with smaller tolerances. In reducing the tolerances all components of RTOL and 
ATOL, and hence of F.WT . should be scaled down uniformly. 

There is no requirement that the same values for ITOL, RTOL, and ATOL be 
used throughout the problem. If during the course of the problem any of these 
parameters is changed, the user should reset ISTATE = 3 before calling LSODE 
again. (ISTATE need not be reset; however, LSODE will not then check the 
legality of the new values.) This option is useful, for example, if the solution 
displays rapid changes in a small subinterval but is relatively smooth elsewhere. 
To accurately track the solution in the rapidly varying region, small values of 
RTOL and ATOL may be required. However, in the smooth regions these 
tolerances could be increased to minimize execution time. 

4.5.8 Selection of Integration and Corrector Iteration Methods (MF) 

The choice of the method flag MF may also require some experimentation. The 
user should consider the nature of the problem and storage requirements. The 
primary consideration regarding MF is stiffness. If the problem is not stiff, the 
best choice is probably MF = 10 (Adams-Moulton (AM) method with functional 
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4. Description of Code Usage 

iteration.) If the problem is stiff to a significant degree, METH should be set 
equal to 2 (table 3.1), and MITER (table 3.2) depends on the structure of the 
Jacobian matrix. If the Jacobian is banded, MITER = 4 (user-supplied analytical 
Jacobian) or 5 (internally generated Jacobian by finite-difference approximations) 
should be used. For either of these two MITER values the user must set values for 
the lower (ML) and upper (MU) half-bandwidths of the Jacobian matrix. The first 
and second elements of the integer work array IWORK must be set equal to ML 
and MU, respectively; that is, IWORK(l) = ML and IWORK(2) = MU. For a full 
matrix MITER should be set equal to 1 (analytical Jacobian) or 2 (internally 
generated Jacobian). If the matrix is significantly diagonally dominant, the choice 
MITER = 3, that is, Jacobi-Newton (JN) iteration using an internally generated 
diagonal approximation for the Jacobian matrix, can be made. To use this 
iteration technique with an analytical Jacobian, set MITER = 4 and ML = MU = 0. 

If the problem is only mildly stiff, the choice METH = 1 (i.e., the AM method) 
may be more efficient than METH = 2 (i.e., the backward differentiation formula 
(BDF) method). For this case experimentation would be necessary to identify the 
optimal METH. If the user has no a priori knowledge regarding the stiffness of 
the problem, one way to determine its nature is to try MF = 10 and examine the 
behavior of both the solution and step size pattern. (It is recommended that some 
upper limit be set for the total number of steps or derivative evaluations to avoid 
excessive run times.) If the typical values of the step size are much smaller than 
the solution behavior would appear to require, for example, more than 100 steps 
are taken over an interval in which the solution changes by less than 1 percent, the 
problem is probably stiff. The degree of stiffness can be estimated from the step 
sizes used and the smoothness of the solution. 

Irrespective of the integration method selected, the least effective iteration 
technique is functional iteration, given by MITER = 0, and the most effective is 
Newton-Raphson (NR), given by MITER = 1 or 2 (4 or 5 for a banded Jacobian 
matrix). Generally JN iteration is somewhere in between. However, storage 
requirements increase in the same order as the effectiveness of the iteration 
technique (see table 4.9), and so trade-off considerations are necessary. For 
reasons of computational efficiency the user is encouraged to provide a routine for 
computing the analytical Jacobian, unless the system is fairly complicated and 
analytical expressions cannot be derived for the matrix elements. The accuracy of 
the Jacobian calculation can be checked by comparison with the J internally 
generated with MITER = 2 or 5. Jacobi-Newton iteration requires considerably 
less storage and execution time per iteration but will be effective only if the 
Jacobian matrix is significantly diagonally dominant. 

The importance of supplying an analytical Jacobian matrix, especially for large 
problems, is illustrated by Radhakrishnan (ref. 37), who studied 12 test problems 
from combustion kinetics. The problems covered a wide range of reaction 
conditions and reaction mechanism size. The effects on solution efficiency of 
(1) METH, (2) the first output station, and (3) optimizing the local error tolerances 
were also examined. 


4.6 Optional Input 

4.5.9 Switching Integration and Corrector Iteration Methods 

The user may specify different values for MF in different subintervals of the 
problem. This option is useful if the problem changes character and is nonstiff in 
some regions and stiff elsewhere. Because stiff problems are usually characterized 
by a nonstiff initial “transient” region, one could use MF = 10 in the initial region 
and then switch to MF = 21 (the BDF method with NR iteration using an 
analytical Jacobian matrix) in the later stiff regime. It is very straightforward to 
change integration methods and corrector iteration techniques. Upon return from 
LSODE the user simply resets MF to the desired new value. The other action 
required is to reset ISTATE = 3 before calling LSODE again. The lengths LRW 
and LIW, respectively, of the arrays RWORK and IWORK depend on MF (see 
tables 4.9 and 4. 10). If different methods are to be used in the course of solving a 
problem, storage corresponding to at least the maximum values of LRW and LIW 
must be allocated. That is, the dimensions of RWORK and IWORK must be set 
equal to at least the largest of the LRW and LIW values, respectively, required by 
the different methods to be used. 


4.6 Optional Input 

In addition to the input parameters whose values are required by the code, the 
user can set values for several other parameters to control both the integration and 
the output from the code. These optional input parameters are given in table 4.6, 
together with their locations and default values. If any of these parameters are 
used, the user must set IOPT = 1 to relay this information to the solver, which will 
examine all optional input parameters and select only those for which nonzero 
values are specified. A value of zero for any parameter will cause its default value 
to be used. Thus to use a subset of the optional inputs, set RWORK(I) = 0.0 and 
IWORK(I) = 0 (I = 5 to 7), and then set parameters of interest to the desired 
(nonzero) values. The variable HO, the step size to be attempted on the first step, 
must indicate the direction of integration. That is, HO must be a positive quantity 
for integration in the forward direction (increasing values of the independent 
variable) and negative otherwise. All other input parameters must be positive 
numbers; otherwise, an error exit will occur. 

To reset any optional input parameter on a subsequent call to LSODE, ISTATE 
must be set equal to 3. IOPT is not altered by LSODE and therefore need not be 
reset. Also because the code does not alter the values in RWORK (5) to RWORK 
(7) and IWORK(5) to IWORK(7), only parameters for which new values are 
required need to be reset. To specify a default value for any parameter for which a 
nondefault value had previously been used, simply load the appropriate location 
in RWORK or IWORK with a zero. Of course, if all variables are to have default 
values, simply reset IOPT = 0. 
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4. Description of Code Usage 

4.6.1 Initial Step Size (HO) 

The sign of the step size HO must agree with the direction of integration; 
otherwise, an error exit will occur. Also, its magnitude should be considerably 
smaller than the average value expected for the problem because the code starts 
the integration with a first-order method. Of course, the integrator tests that the 
given step size does produce a solution that satisfies the local error test and, if 
necessary, decreases it (in magnitude). The only test made on the magnitude of 
HO prior to taking the first step is that it does not exceed the user-supplied value 
for HMAX, the maximum absolute step size allowed for the problem. 

4.6.2 Maximum Step Size (HMAX) 

The user may have to specify a finite value for HMAX (default value, if the 
solution is characterized by rapidly varying transients between long smooth 
regions. If the step size is too large, the solver may skip over the fine detail that 
the user may be (primarily) interested in. An example of this behavior is the 
buildup of ozone and oxygen atom concentrations in the presence of sunlight 
(ref. 17). 

4.6.3 Maximum Method Order (MAXORD) 

The optional input parameter MAXORD, the maximum method order to be 
attempted on any step, should not exceed the default value — 12 for the AM 
method and 5 for the BDF method. If it does, it will be reduced to the default 
value. Also, in the course of solving the problem, if MAXORD is decreased to a 
value less than the current method order, the latter quantity will be reduced to the 
new MAXORD. 

The maximum method order has to be restricted to a value less than the default 
value for stiff problems when the eigenvalues of the Jacobian matrix are close to 
the imaginary axis; that is, the solution is highly oscillatory. In such a situation 
the BDF method of high order (> 3) has poor stability characteristics and, as the 
stability plots in Gear (ref. 10) show, the unstable region grows as the order is 
increased. For this reason MAXORD should be set equal to 3 unless the 
eigenvalues are imaginary; that is, Re(?t;) = 0 and Im (X/) * 0, where Re(A,,) and 
lm(Xj ) are the real and imaginary parts of the rth eigenvalue. In this case the 
value MAXORD = 2 should be used. 


4.7 Optional Output 

The user is usually primarily interested in the numerical solution and the 
corresponding value of the independent variable. These quantities are always 
returned in the call variables Y and T. In addition, several optional output 
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4.8. Other Routines 


quantities that contain information about the integration are returned by LSODE. 
These quantitites are given in tables 4.7 and 4.8, together with their locations. 
Some of these quantities give a measure of the computational work required and 
may, for example, help the user decide if the problem is stiff or if the right method 
is being used. Other output quantities will, in the event of an error exit, help the 
user either set legal values for some parameters or identify the reason for repeated 
convergence failures or local error test failures. 


4.8 Other Routines 

To gain additional capabilities, the user can access the following subroutines 

included in the LSODE package: INTDY, SRCOM, XSETF, and XSETUN. 

Among these, only INTDY is used by LSODE. 

4,8.1 Interpolation Routine (Subroutine INTDY) 

The subroutine INTDY provides derivatives of Y, up to the current order, at a 

specified point T and may be called only after a successful return from LSODE. 

The call to this routine takes the form 

CALL INTDY (T, K, RWORK(2 1 ), NYH, DKY, IFL AG) . 

where T, K, RWORK(21), and NYH are input parameters and DKY and IFLAG 

are output parameters. The arguments to INTDY are defined as follows: 

T Value of independent variable at which the results are required. 

For the results to be valid T must lie in the interval [(TCUR - 
HU),TCUR], where TCUR and HU are defined in table 4.7. 

K Integer that specifies the desired derivative order and must satisfy 

0 < K < current method order NQCUR (see table 4.7 for location 
of this quantity). Now, because the method order is never less 
than 1, the first derivative dY/d^ can always be obtained by 
calling INTDY. 

RWORK(2 1 ) Base address of the Nordsieck history array (see table 4.8). 

NYH Number of ODE’s used on the first call to LSODE. If the number 

of ODE’s is decreased during the course of the problem, NYH 
should be saved. An alternative way of obtaining NYH is to 
include the common block LS0001 in the subprogram calling 
INTDY. LSODE saves NYH in LS0001 as the 232nd word — the 
14th integer word after 218 real words (see table 3.6). 
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DKY Array of length N that contains the Kth derivative of Y at T. The 

subprogram calling INTDY must include a DIMENSION statement 
for DKY if NYH > 1 . Alternatively, to save storage, DKY can be 
replaced with RWORK(LSAVF)— see section 4.3. 

IFLAG An error flag with following values and meanings: 

0 Both T and K were legal. 

-1 Illegal value was specified for K. 

-2 Illegal value was specified for T. 

4.8.2 Using Restart Capability (Subroutine SRCOM) 

The subroutine SRCOM is useful if one is either alternating between two or 
more problems being solved by LSODE or interested in interrupting a run and 
restarting it later. The latter situation may arise, for example, if one is interested in 
steady-state values with no a priori knowledge of the required integration interval. 
The run may be stopped periodically, the results examined and, if necessary, the 
integration continued. This procedure is clearly more economical than making 
repeated runs on the same problem with, say, increasing values of TOUT. To 
exploit the capability of stopping and then continuing the integration, the user 
must save and then restore the contents of the common blocks LS0001 and 
EH0001. This information can be stored and restored by calling SRCOM. The 
call to this routine takes the form 

CALL SRCOM (RSAV, ISAV, JOB) 

where RSAV must be declared as a real array of length 218 or more in the calling 
subprogram and ISAV as an integer array of length 41 or more and JOB is an 
integer flag whose value (= 1 or 2) indicates the action to be performed by 
SRCOM as follows: JOB = 1 means “save the contents of the two common 
blocks,” and JOB = 2 means “restore this information.” 

Thus to store the contents of EH0001 and LS0001, SRCOM should be called as 
follows; 


CALL SRCOM (RSAV, ISAV, 1) 

Upon return from SRCOM, RSAV and ISAV will contain, respectively, the 218 
real and 39 integer words that together make up the common block LS0001 . The 
40th and 41st elements of ISAV will contain the two integer words MESFLG and 
LUNIT in the common block EH0001 (table 3.6). The lengths and contents of the 
arrays RWORK and IWORK must also be saved. The lengths LENRW and 
LENIW required for the arrays RWORK and IWORK are saved by LSODE as the 
17th and 18th elements, respectively, of the array IWORK (see table 4.7). 

To continue the integration, the arrays RWORK and IWORK and the contents 
of the common blocks LS0001 and EH0001 must be restored. The common block 
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contents are restored by using the previously saved arrays RSAV and ISAV and 
calling the routine SRCOM as follows: 

CALL SRCOM (RSAV, ISAV, 2) 

The user should then set values for the input parameters required by LSODE, and 
the integration can be continued by calling this routine. Note, in particular, that 
ISTATE must be set equal to 2 or 3 to inform LSODE that the present call is a 
continuation one for the problem (see table 4.3). 

4.8.3 Error Message Control (Subroutines XSETF and XSETUN) 

To reset the value of the logical unit number LUNIT for output of messages 
from the code, the routine XSETUN should be called as follows: 

CALL XSETUN (LUN) 

where LUN is the new value for LUNIT. Action is taken only if the specified 
value is greater than zero. 

The value of the flag MESFLG, which controls whether messages from the 
code are printed or not, may be reset by calling subroutine XSETF as follows: 

CALL XSETF (MFLAG) 

where MFLAG is the new value for MESFLG. The legal values for MFLAG are 
0 and 1. Specifying any other value will result in no change to the current value 
of MESFLG. Setting MFLAG = 0 does carry the risk of losing valuable information 
through error messages from the integrator. 


4.9 Optionally Replaceable Routines 

If none of the error control options included in the code are suitable, more 
general error controls can be obtained by substituting user-supplied versions of 
the routines EWSET and/or VNORM (table 3.3). Both routines are concerned 
with measuring the local error. Hence any replacement may have a major impact 
on the performance of the code. We therefore recommend that modifications be 
made only if absolutely necessary, and that too with great caution. Also the effect 
of the changes and the accuracy of the programming should be studied on some 
simple problems. 

4.9.1 Setting Error Weights (Subroutine EWSET) 

The subroutine EWSET sets the array of error weights EWT, equation (4.1). 
This routine takes the form 
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SUBROUTINE EWSET (N, ITOL, RTOL, ATOL, YH, EWT) 

where N is the current value of the number of ODE’s; ITOL, RTOL, ATOL, and 
EWT have been defined previously; and YH contains the current Nordsieck 
history array, that is, the current solution vector YCUR and its NQ scaled 
derivatives, where NQ is the current method order. On the first call to EWSET 
from the routine LSODE, YCUR is the same as the Y array (which then contains 
the initial values supplied by the user); thereafter the two arrays may be different. 

The error weights {EWT/} are used in the local truncation error test, which 
requires that the rms norm of dj/EWT, be 1 or less. Here, di is the estimated local 
error in 7/. The above norm is computed in the routine VNORM (discussed in 
section 4.9.2) to which the EWT array is passed. 

If the user replaces the current version of EWSET, the new version must return 
in each EWT/ 0’ = 1,...,N) a positive quantity for comparison with d;. This routine 
is called by the routine LSODE only (tables 3.4 and 3.5). However, in addition to 
its use in the local truncation error test (which is performed in the routine 
STODE), EWT is used (1) by the routine LSODE in computing the initial step 
size HO and the optional output integer IMXER (table 4.7) and (2) by the routine 
PREPJ in computing the increments in solution vector for the difference quotient 
Jacobian matrix (MITER = 2 or 5, table 3.2) and for the diagonal approximation 
to the Jacobian matrix (MITER = 3). The base address for EWT in the array 
RWORK is LEWT, which is the 222nd word (the 4th integer word after 218 real 
words) in the common block LS0001 . 

If the user’s version of EWSET uses current values of the derivatives of i, they 
can be obtained from YH, as described later. Indeed, derivatives of any order, up 
to NQ, can be found from YH, whose base address in RWORK is LYH (= 21), the 
221st word (the 3rd integer word past 218 real words) in LS0001. The array YH is 
of length NYH(NQ + 1), where NYH is the value of N on the first call to LSODE. 
The first N elements correspond exactly to the YCUR array. The remaining terms 
contain scaled derivatives of YCUR. For example, the N elements J*NYH + 1 to 
J*NYH + N (J = 0,1,. ..,NQ) contain the Jth scaled derivative H J Y^/J!, where H is 
the current value of the step size. On the first call to EWSET, before any 
integration is done, H is (temporarily) set equal to 1 .0. Thereafter its value may be 
determined from LS0001, where it is the 212th real word. This common block 
also contains NYH as the 232nd word (the 14th integer word past 218 real words) 
and NQ as the 253rd word (the 35th integer word past 218 real words). Thus if the 
user wishes to use the Jth derivative in EWSET, it may be obtained by including 
the following statements: 

SUBROUTINE EWSET (N, ..., YH, ..., EWT) 

REAL (or DOUBLE PRECISION) YH, EWT, RLS, H, ... 
INTEGER N, JLS, NQ, NYH, ... 

DIMENSION YH(1), EWT(l), ... in FORTRAN 66 
DIMENSION YH(*), EWT(*), ... in FORTRAN 77 
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COMMON/LS0001/RLS(218), ILS(39) 

NQ = ILS(35) 

NYH = ILS(14) 

H = RLS(212) 

The Jth derivative (0 < J < NQ) is then given by 


Y ( J ) - 
M ” 


J! 


yh(j*nyh+i) 

H 7 


1=1 N. 


(4.2) 


where y| J) is the Jth derivative of Yj. The routine must include a data type 
declaration and a DIMENSION statement for To save on storage, these 
values may be stored temporarily in the vector EWT. 

4.9.2 Vector-Norm Computation (Function VNORM) 

The real (or double precision) function routine VNORM computes the weighted 
root-mean-square (rms) norm of a vector. It is used as follows: 

D = VNORM (N, V, W) 

where N is the length of the real arrays V, which contains the vector, and W, which 
contains the weights. Upon return from VNORM, D contains the weighted rms- 
norm 



n y v . v 


,=,i w «v 


This routine is used by STODE to compute the weighted rms norm of the 
estimated local error. STODE also uses information returned by VNORM to 
perform the corrector convergence test and to compute factors that determine if 
the method order should be changed. Other routines that access VNORM are 
LSODE, to compute the initial step size HO, and PREPJ, to compute the increments 
in the solution vector for generating difference quotient Jacobians (MITER = 2 or 
5, table 3.2). 

If the user replaces the routine VNORM, the new version must return a positive 
quantity in VNORM, suitable for use in local error and convergence testing. The 
weight array W can be used as needed, but it must not be altered in VNORM. For 
example, the max-norm, that is, maxIV/W,!, satisfies this requirement, as does a 
norm that ignores some components of V. The latter procedure has the effect of 
suppressing error control on the corresponding components ofY. 
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4.10 Overlay Situation 

If LSODE is to be used in an overlay situation, the user must declare the 
variables in the call sequence to LSODE and in the two internal common blocks 
LS0001 and EH0001 in the MAIN program to ensure that their contents are 
preserved. The common block LS0001 is of length 257 (218 real or double- 
precision words followed by 39 integer words), and EH0001 contains two integer 
words (see table 3.6). 


4.11 Troubleshooting 

In this section we present a brief discussion of the corrective actions that may 
be taken in case of difficulty with the code. If the execution is terminated 
prematurely, the user should examine the error message and the value of ISTATE 
returned by LSODE (table 4.4). We therefore recommend that the current value of 
MESFLG not be changed, at least until the user has gained some experience with 
the code. The legality of every input parameter, both required and optional, is 
checked. If illegal input is detected by the code, it returns to the calling subprogram 
with ISTATE = -3. The error message will be detailed and will make clear what 
corrective actions to take. If the illegal input is caused by a request for too much 
accuracy, the user should examine the value of TOLSF returned in RWORK(13) 
(table 4.7) and make necessary adjustments to RTOL and ATOL, as described in 
section 4.5.7. If an excessive accuracy requirement is detected during the course 
of solving the problem, the value ISTATE = —2 is returned. To continue the 
integration, make the adjustments mentioned above, set ISTATE = 3, and call 
LSODE again. 

Another difficulty related to accuracy control may be encountered if pure 
relative error control for, say, the ith variable is specified (i.e., ATOL, = 0). If this 
solution component vanishes, the error test cannot be applied. In this situation the 
value ISTATE = -6 is returned to the calling subprogram. The error message 
identifies the component causing the difficulty. To continue integrating, reset 
ATOL for this component to a nonzero value, set ISTATE = 3, and call LSODE 
again. 

If more than MXSTEP (default value, 500) integration steps are taken on a 
single call to LSODE without completing the task, the error return ISTATE = -1 is 
made. The problem might be the use of an inappropriate integration method or 
iteration technique. The use of ME = 10 (or 20) on a stiff problem is one example. 
The user should, as described previously under the selection of MF (section 
4.5.8), verify that the value of MF is right for the problem. Very stringent accuracy 
requirements may also cause this difficulty. Another possibility is that pure 
relative error control has been specified but most, or all, of the |yj| are very small 
but nonzero. Finally, the solution may be varying very rapidly, forcing the 
integrator to select very small step sizes, or the integration interval may be very 
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long relative to the average step size. To continue the integration, simply reset 
ISTATE = 2 and call LSODE again — the excess step counter will be reset to zero. 
To prevent a recurrence of the error, the value of MXSTEP can be increased, as 
described in section 4.6. If this action is taken between calls to LSODE, ISTATE 
must be set equal to 3 before LSODE is called again. Irrespective of when 
MXSTEP is increased, IOPT should be set equal to 1 before the next call to 
LSODE. 

If the integrator encounters either repeated local error test failures or any local 
error test failure with a step size equal to the user-supplied minimum value HMIN 
(table 4.6), LSODE returns with ISTATE = -A. The difficulty could be caused by 
a singularity in the problem or by inappropriate input. The user should check 
subroutines F and JAC for errors. If none is found, it may be necessary to monitor 
intermediate quantities. The component IMXER causing the error test failure is 
returned as IWORK(16) (table 4.7). The values Y(IMXER), RTOL(IMXER), 
ATOL(IMXER), and ACOR(IMXER) (see table 4.8) should be examined. If pure 
relative error control had been specified for this component, very small but 
nonzero values of Y(IMXER) may cause the difficulty. 

These checks should also be made if the integration fails because of either 
repeated corrector convergence test failures or any such failure with a step size 
equal to HMIN. In this case LSODE returns the value ISTATE = -5 along with a 
value for IMXER defined above. If an analytical Jacobian is being used, it should 
be checked for errors. The accuracy of the calculation can also be checked by 
comparing J with that generated internally. Another reason for this failure may be 
the use of an inappropriate MITER, for example, MITER = 3 for a problem that 
does not have a diagonally dominant Jacobian. It may be helpful to try different 
values for MITER and monitor the successive corrector estimates stored as the Y 
array in subroutine STODE. 

In addition to the error messages just discussed, a warning message is printed if 
the step size H becomes so small that T + H = T on the computer, where T is the 
current value of the independent variable. This error is not considered fatal, and 
so the execution is not terminated nor is a return made to the calling subprogram. 
No action is required by the user. The warning message is printed a maximum 
number of MXHNIL (default value, 10) times per problem. The user can change 
the number of times the message is printed by resetting MXHNIL, as discussed in 
section 4.6. To indicate the change to LSODE, the parameter IOPT must be set 
equal to 1 before LSODE is called. 
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Chapter 5 

Example Problem 

5.1 Description of Problem 

In this chapter we demonstrate the use of the code by means of a simple stiff 
problem taken from chemical kinetics. The test case, described elsewhere (refs. 
17, 28, and 38), consists of three chemical species participating in three irreversible 
chemical reactions at constant density and constant temperature: 

Sj -> S 2 , (5.1) 


*2 

s 2 +fi 3 -> fij +s 3 , (5.2) 

*3 

S 2 +S 2 2 3 +S 3 , (5.3) 

with k\ = 4xl0 -2 , ki = 10 4 , and k$ = 1.5xl0 7 . In reactions (5.1) to (5.3), S* is the 
chemical symbol for the ith species, the arrows denote the directions of the 
reactions (the single arrow for each reaction means that it takes place in the 
indicated direction only), and the [kj] are the specific rate coefficients for the 
reactions. The units of kj depend on reaction type (e.g., ref. 39). If y ; denotes the 
molar concentration of species i, that is, moles of species i per unit volume of 
mixture, the governing ODE’s are given by 



-0.04 y, H-IO 4 }^, 


(5.4) 


dy 

dt 


h = 0.04 -10 4 y 2 3' 3 -3xl0 7 y 2 >' 2 , 


(5.5) 
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^3 

dt 


= 3x10 7 >'2) , 2’ 


where t is time in seconds. The initial conditions are 


(5-6) 


y x (t = 0) = 1 ; y 2 (t = 0) = y 3 (t = 0) = 0. (5.7) 

The example problem is interesting because the reaction rate coefficients vary 
over nine orders of magnitude. Also it can be quite easily verified that at steady 
state, that is, as t -> 0, >2 -» 0, and J 3 -> 1. To study the evolution of the 

chemical system, including the approach to the final state, we integrate the ODE’s 
up to t = 4xl0 10 s, generating output at t = 0.4x10” s (n = 0,1,. ..,11). 

5.2 Coding Required To Use LSODE 

5.2.1 General 

All of the coding required to solve the example problem with LSODE is 
included (in the form of comment statements) in the package supplied to the user. 
The MAIN program that calls LSODE and manages output is given in figure 5. 1 . 
Figure 5.2 lists the subroutine that computes the derivatives. Because a value of 
MITER = 1 is used (fig. 5.1), a routine that computes the analytical Jacobian 
matrix is required. This routine is given in figure 5.3. The names used for the 
derivative and Jacobian matrix subroutines are, respectively, FEX and JEX. 
Therefore these names are used as arguments in the call to LSODE and declared 
EXTERNAL in the MAIN program (fig. 5.1). 

5.2.2 Selection of Parameters 

Because the problem is stiff, the choice METH = 2 is made. For the same 
reason functional iteration, that is, MITER = 0, is rejected. It is straightforward to 
compute the analytical Jacobian matrix, which should be used for reasons of 
efficiency. In any case, the choice MITER = 3, that is, Jacobi-Newton iteration, 
must not be made because the Jacobian matrix is not diagonally dominant. The 
choice MITER = 4 with ML = 1 and MU = 2 could be made but will require more 
storage than MITER = 1 (see table 4.9). More importantly the computational 
overhead for the LU-decomposition of the iteration matrix is more for MITER = 4 
than for MITER = 1 . Hence the value MF = 21 is used. 

The number NEQ of ODE’s is equal to the number (= 3) of chemical species. 
To minimize storage, the lengths LRW and LIW of the work arrays RWORK and 
IWORK are set equal to their minimum required values. According to the 
formulas given in tables 4.9 and 4.10 for MF = 21, these lengths are as follows: 
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EXTERNAL FEX, JEX 

DOUBLE PRECISION ATOL, RWORK, RTOL, T, TOUT, Y 
DIMENSION Y(3), AT0L(3), RW0RK(58), IW0RK{23) 

NEQ • 3 
Y(l) » 1.D0 
Y(Z) - 0.D0 
Y(3) - 0.D0 
T • 0.D0 
TOUT * .4DO 
ITOL • 2 
RTOL * l.D-4 
ATOL(l) * l.D-6 
AT0L(2) ■ l.D-10 
AT0L(3) « l.D-6 
ITASK ■ 1 
ISTATE - 1 
IOPT - 0 
LRW ■ 58 
LIW • 23 
MF » 21 

DO 40 I OUT - 1,12 

CALL LSODE (FEX, NEQ, Y,T, TOUT, ITOL, RTOL, ATOL, ITASK, ISTATE, 

1 IOPT, RWORK, LRW, IWORK, LIW, JEX, MF) 

WRITE (3 , 20) T , Y (1) , Y (2) , Y(3) 

20 F0RMAT(7H AT T -,E12.4,6H Y -.3E15.7) 

IF (ISTATE .LT. 0) 60 TO 80 
40 TOUT ■ TOUTMO.DO 

WRITE(3,60)IW0RK(ll),IWORK(12) ,IW0RK(13) 

60 FORMAT (/12H NO. STEPS -,I4,11H NO. F-S -,I4,11H NO. J-S «,I4) 
STOP 

80 WRITE (3, 90) 1ST ATE 

90 FORMAT(///22H ERROR HALT.. ISTATE »,I3) 

STOP 

ENO 


Figure 5.1 .—Listing of MAIN program for example problem. 


SUBROUTINE FEX (NEQ, T, Y, YDOT) 

DOUBLE PRECISION T, Y, YDOT 
DIMENSION Y (3) , YD0T(3) 

YDOT(l) - -,04D0*Y(1) + 1 . D4*Y (2) *Y (3) 
YDOT 3 « 3.D7*Y(2)*Y(2) 

YDOT (2) - -YDOT(l) - YD0T(3) 

RETURN 

END 

Figure 5.2. — Listing of subroutine (FEX) that 
computes derivatives for example problem. 


LRW = 22 + 3(5 + 1) + 3(3) + 3 2 = 58 
and 


LIW = 20 + 3 = 23. 

Selection of the error tolerances requires some explanation. A scalar RTOL is 
used because the same number of significant figures is acceptable for all 
components. However, because yi is expected to be much smaller than both y\ 
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SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD) 
DOUBLE PRECISION PD, T, Y 
DIMENSION Y{3), PD(NRPD,3) 


PD 

PD 

PO 

PD 

PD 

PD 

PO 


1.1) 

1 . 2 ) 

1.3 

2,1 

2.3 
3,2 
2 , 2 ) 


RETURN 

END 


-.0400 
l.D4*Y( 
l.D4*Y( 
.0400 
-PD(1,3) 
6.D7*Y(2) 
-P0(1.2) - 


PD(3,2) 


Figure 5.3.— Listing of subroutine (JEX) that computes 
analytical Jacobian matrix for example problem. 


and an array ATOL, with ATOL(2) much smaller than both ATOL(l) and 
ATOL(3), is used. For these choices of the RTOL and ATOL types, table 4.2 gives 
ITOL = 2. Pure relative error control cannot be used because the initial values of 
both >’2 and are zero and, as t — > y\ 0 and j2 ~ * 0- Pure absolute error 
control should not be used because of the widely varying orders of magnitude of 
the {y/}. Note that because a scalar RTOL is used, the MAIN program does not 
require a DIMENSION statement for this variable. 

The remainder of the program calling LSODE is straightforward and self- 
explanatory. Because the output value for ISTATE is equal to 2 for a normal 
return from LSODE and no parameter (except TOUT) is reset between calls to 
LSODE, ISTATE does not have to be reset. 


5.3 Computed Results 


The output from the program, obtained on the Lawrence Livermore Laboratory’s 
CDC-7600 computer using single-precision arithmetic, is given in figure 5,4. In 
addition to the results at the specified times, values for the following parameters, 
which give a measure of the computational work required to solve the problem, 
are printed at the end: total number of integration steps (STEPS), total number of 
derivative evaluations (F— S), and total number of Jacobian matrix evaluations and 
LU-decompositions of the iteration matrix (J-S). 


AT T 

a 

4.0000E-01 

Y 

- 9.851726E-01 

3.386406E-05 

1.479357E-02 

AT T 

■ 

4.0000E+00 

Y 

- 9.055142E-01 

2.240418E-05 

9.446344E-02 

AT T 

■ 

4.0000E+01 

Y 

- 7.158050E-0! 

9.184616E-06 

2.841858E-01 

AT T 

m 

4.0000E+02 

Y 

* 4.504846E-01 

3.222434E-06 

5.495122E-01 

AT T 

m 

4.0000E+03 

Y 

- 1.831701E-01 

8.940379E-07 

8. 168290E-01 

AT T 

M 

4.0000E+04 

Y 

- 3.897016E-02 

1.621193E-07 

9.610297E-01 

AT T 

u 

4.0000E+05 

Y 

- 4.935213E-03 

1.983756E-08 

9.950648E-01 

AT T 

M 

4.0000E+06 

Y 

- 5.159269E-04 

2.064759E-09 

9.994841E-01 

AT T 

m 

4.0000E+07 

Y 

- 5. 306413E-05 

2.122677E-10 

9.999469E-01 

AT T 

m 

4.0000E+08 

Y 

» 5.494529E-D6 

2.197824E-11 

9.999945E-01 

AT T 

a 

4.0000E+09 

Y 

- 5.129458E-07 

2.051784E-12 

9.999995E-01 

AT T 

m 

4.0000E+10 

Y 

- -7.170592E-08 ■ 

-2.868236E-13 

1.000 00 OE+OO 

NO. STEPS » 330 NO. 

F-S * 405 NO. J-S 

- 69 



Figure 5.4. — Output from program for example problem. 
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Code Availability 

The present version of LSODE, dated March 30, 1987, is available in single or 
double precision. The code has been successfully executed on the following 
computer systems: Lawrence Livermore Laboratory’s CDC-7600, Cray-1, and 
Cray-X/MP; NASA Lewis Research Center’s IBM 370/3033 using the TSS 
operating sytem (OS), Amdahl 5870 using the VM/CMS OS and the UTS OS, 
Cray-X/MP/2/4 using the COS and UNICOS operating sytems and the CFT and 
CFT77 compilers, Cray-Y/MP/8/6128 using UNICOS 6.0 and CFT77, Alliant 
FX/S, Convex C220 minicomputer using the Convex 8.0 OS, and VAX 
1 1/750, 1 1/780, 1 1/785, 6320, 6520, 8650, 8800, and 9410; NASA Ames Research 
Center’s Cray-2 and Cray-Y/MP using the UNICOS operating system and the 
CFT77 compiler; the Sun SPARCstation 1 using the Sun 4.0 OS; the IBM RISC 
System/6000 using the AIX 3. 1 OS and the XLF and F77 compilers; several IRIS 
workstations using the IRIX 4.0.1 OS and F77 compiler; and various personal 
computers under various systems. 

The LSODE package is one of Five solvers included in the ODEPACK collection 
of software for ordinary differential equations (ref. 2). The official distribution 
center for ODEPACK is the Energy Science and Technology Software Center at 
Oak Ridge, Tennessee. (ESTSC supersedes NESC, the National Energy Software 
Center at Argonne National Laboratory, in this activity.) Both single- and double- 
precision versions of the collection are available. Additional details regarding 
code availability and procurement can be obtained from 

Energy Science and Technology Software Center 

P.O. Box 1020 

Oak Ridge, TN 37831-1020 

Telephone: (615) 576-2606 

The ODEPACK solvers can also be obtained through electronic mail by accessing 
the NETLIB collection of mathematical software (ref. 40). Both single- and 
double-precision versions of ODEPACK are contained in NETLIB. Detailed 
instructions on how to access and use NETLIB are given by Dongarra and Grosse 
(ref. 40). 
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